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ESSENTIALLY NON-OSCILLATORY AND WEIGHTED ESSENTIALLY 
NON-OSCILLATORY SCHEMES FOR HYPERBOLIC CONSERVATION LAWS 


CHI- WANG SHU * 

Abstract. In these lecture notes we describe the construction, analysis, and application of ENO (Es- 
sentially Non-Oscillatory) and WENO (Weighted Essentially Non- Oscillatory) schemes for hyperbolic con- 
servation laws and related Hamilton- Jacobi equations. ENO and WENO schemes are high order accurate 
finite difference schemes designed for problems with piecewise smooth solutions containing discontinuities. 
The key idea lies at the approximation level, where a nonlinear adaptive procedure is used to automatically 
choose the locally smoothest stencil, hence avoiding crossing discontinuities in the interpolation procedure as 
much as possible. ENO and WENO schemes have been quite successful in applications, especially for prob- 
lems containing both shocks and complicated smooth solution structures, such as compressible turbulence 
simulations and aeroacoustics. 

These lecture notes are basically self-contained. It is our hope that with these notes and with the help of 
the quoted references, the reader can understand the algorithms and code them up for applications. Sample 
codes are also available from the author. 

Key words, essentially non- oscillatory, conservation laws, high order accuracy 

Subject classification. Applied and Numerical Mathematics 

1. Introduction, ENO (Essentially Non-Oscillatory) schemes started with the classic paper of Harten, 
Engquist, Osher and Chakravarthy in 1987 [38], This paper has been cited at least 144 times by early 1997, 
according to the ISI database. The Journal of Computational Physics decided to republish this classic paper 
as part of the celebration of the journal’s 30th birthday [68], 

Finite difference and related finite volume schemes are based on interpolations of discrete data using 
polynomials or other simple functions. In the approximation theory, it is well known that the wider the 
stencil, the higher the order of accuracy of the interpolation, provided the function being interpolated is 
smooth inside the stencil. Traditional finite difference methods are based on fixed stencil interpolations. For 
example, to obtain an interpolation for cell i to third order accuracy, the information of the three cells i — 1, 
i and i + 1 can be used to build a second order interpolation polynomial. In other words, one always looks 
one cell to the left, one cell to the right, plus the center cell itself, regardless of where in the domain one 
is situated. This works well for globally smooth problems. The resulting scheme is linear for linear PDEs, 
hence stability can be easily analyzed by Fourier transforms (for the uniform grid case). However, fixed 
stencil interpolation of second or higher order accuracy is necessarily oscillatory near a discontinuity, see 
Fig. 2.1, left, in Sect. 2.2. Such oscillations, which are called the Gibbs phenomena in spectral methods, do 
not decay in magnitude when the mesh is refined. It is a nuisance to say the least for practical calculations, 
and often leads to numerical instabilities in nonlinear problems containing discontinuities. 

Before 1987, there were mainly two common ways to eliminate or reduce such spurious oscillations near 
discontinuities. One way was to add an artificial viscosity. This could be tuned so that it was large enough 

* Division of Applied Mathematics, Brown University, Providence, RI 02912 (e-mail: shu@cfm.brown.edu). Research of 
the author was partially supported by NSF grants DMS-9500814, ECS- 9214488, ECS- 9627849 and INT- 9601 084, ARO grants 
DAAH04-94-G-0205 and DAAG55-97-1-0318, NASA Langley grant NAG-1-1145 and Contract NAS1-19480 while in residence 
at ICASE, NASA Langley Research Center, Hampton, VA 23681-0001, and AFOSR grant F49620-96- 1-0150. 
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near the discontinuity to suppress, or at least to reduce the oscillations, but was small elsewhere to maintain 
high-order accuracy. One disadvantage of this approach is that fine tuning of the parameter controlling 
the size of the artificial viscosity is problem dependent. Another way was to apply limiters to eliminate 
the oscillations. In effect, one reduced the order of accuracy of the interpolation near the discontinuity 
(e.g. reducing the slope of a linear interpolant, or using a linear rather than a quadratic interpolant near 
the shock). By carefully designing such limiters, the TVD (total variation diminishing) property could be 
achieved for nonlinear scalar one dimensional problems. One disadvantage of this approach is that accuracy 
necessarily degenerates to first order near smooth extrema. We will not discuss the method of adding explicit 
artificial viscosity or the TVD method in these lecture notes. We refer to the books by Sod [75] and by 
LeVeque [52], and the references listed therein, for details. 

The ENO idea proposed in [38] seems to be the first successful attempt to obtain a self similar (i.e. no 
mesh size dependent parameter), uniformly high order accurate, yet essentially non-oscillatory interpolation 
(i.e. the magnitude of the oscillations decays as 0(Ax fc ) where k is the order of accuracy) for piecewise 
smooth functions. The generic solution for hyperbolic conservation laws is in the class of piecewise smooth 
functions. The reconstruction in [38] is a natural extension of an earlier second order version of Harten and 
Osher [37]. In [38], Harten, Engquist, Osher and Chakravarthy investigated different ways of measuring local 
smoothness to determine the local stencil, and developed a hierarchy that begins with one or two cells, then 
adds one cell at a time to the stencil from the two candidates on the left and right, based on the size of 
the two relevant Newton divided differences. Although there are other reasonable strategies to choose the 
stencil based on local smoothness, such as comparing the magnitudes of the highest degree divided differences 
among all candidate stencils and picking the one with the least absolute value, experience seems to show 
that the hierarchy proposed in [38] is the most robust for a wide range of grid sizes, Ax, both before and 
inside the asymptotic regime. 

As one can see from the numerical examples in [38] and in later papers, many of which being mentioned 
in these lecture notes or in the references listed, ENO schemes are indeed uniformly high order accurate and 
resolve shocks with sharp and monotone (to the eye) transitions. ENO schemes are especially suitable for 
problems containing both shocks and complicated smooth flow structures, such as those occurring in shock 
interactions with a turbulent flow and shock interaction with vortices. 

Since the publication of the original paper of Harten, Engquist, Osher and Chakravarthy [38], the 
original authors and many other researchers have followed the pioneer work, improving the methodology and 
expanding the area of its applications. ENO schemes based on point values and TVD Runge-Kutta time 
discretizations, which can save computational costs significantly for multi space dimensions, were developed in 
[69] and [70]. Later biasing in the stencil choosing process to enhance stability and accuracy were developed 
in [28] and [67]. Weighted ENO (WENO) schemes were developed, using a convex combination of all 
candidate stencils instead of just one as in the original ENO, [53], [43]. ENO schemes based on other than 
polynomial building blocks were constructed in [40], [16]. Sub-cell resolution and artificial compression to 
sharpen contact discontinuities were studied in [35], [83], [70] and [43]. Multidimensional ENO schemes 
based on general triangulation were developed in [1]. ENO and WENO schemes for Hamilton- Jacobi type 
equations were designed and applied in [59], [60], [50] and [45]. ENO schemes using one-sided Jocobians for 
field by field decomposition, which improves the robustness for calculations of systems, were discussed in 
[25]. Combination of ENO with multiresolution ideas was pursued in [7]. Combination of ENO with spectral 
method using a domain decomposition approach was carried out in [8]. On the application side, ENO and 
WENO have been successfully used to simulate shock turbulence interactions [70], [71], [2]; to the direct 
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simulation of compressible turbulence [71], [80], [49]; to relativistic hydrodynamics equations [24]; to shock 
vortex interactions and other gas dynamics problems [12], [27], [43]; to incompressible flow problems [26], 
[31]; to viscoelasticity equations with fading memory [72]; to semi-conductor device simulation [28], [41], [42]; 
to image processing [59], [64], [73]; etc. This list is definitely incomplete and may be biased by the author’s 
own research experience, but one can already see that ENO and WENO have been applied quite extensively 
in many different fields. Most of the problems solved by ENO and WENO schemes are of the type in which 
solutions contain both strong shocks and rich smooth region structures. Lower order methods usually have 
difficulties for such problems and it is thus attractive and efficient to use high order stable methods such as 
ENO and WENO to handle them. 

Today the study and application of ENO and WENO schemes are still very active. We expect the 
schemes and the basic methodology to be developed further and to become even more successful in the 
future. 

In these lecture notes we describe the construction, analysis, and application of ENO and WENO schemes 
for hyperbolic conservation laws and related Hamilton- Jacobi equations. They are basically self-contained. 
Our hope is that with these notes and with the help of the quoted references, the readers can understand 
the algorithms and code them up for applications. Sample codes are also available from the author. 

2. One Space Dimension. 

2.1. Reconstruction and Approximation in ID. In this section we concentrate on the problems 
of interpolation and approximation in one space dimension. 

Given a grid 

(2.1) a = x i < xi < ... < x N _ i < = by 

We define cells, cell centers, and cell sizes by 

h = yXi — — i 

(2.2) Axi = - x { _iy i = 1,2,..., AT. 

We denote the maximum cell size by 

(2.3) Ax = max Ax^ . 

l<i<N 

2.1.1. Reconstruction from cell averages. The first approximation problem we will face, in solving 
hyperbolic conservation laws using cell averages (finite volume schemes, see Sect. 2.3.1), is the following 
reconstruction problem [38]. 

Problem 2.1. One dimensional reconstruction. 

Given the cell averages of a function v{x)\ 

(2.4) v i = i — 1,2,..., AT , 

find a polynomial p*(x), of degree at most k — 1, for each cell Ji, such that it is a fc-th order accurate 
approximation to the function v(x) inside /*: 

(2.5) pi(x) = v(x) -f 0(Ax k ), x E Ii, i = l,...,iV. 
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In particular, this gives approximations to the function v{x) at the cell boundaries 
(2.6) v7 + i=Pi(x i+ 1 ), v+_^ = Pi (x i _i), i = , 

which are fc-th order accurate: 

(2-7) v~ + ^ =v(x i+ i) + 0(Ax fc ), vf_ i =v(ar i _i) + 0(Ax fc ), i = l,.:.,N. 

The polynomial Pi(x) in Problem 2.1 can be replaced by other simple functions, such as trigonometric 
polynomials. See Sect. 4.1.3. 

We will not discuss boundary conditions in this section. We thus assume that Vi is also available for 
i < 0 and i > TV if needed. 

In the following we describe a procedure to solve Problem 2.1. 

Given the location I t and the order of accuracy fc, we first choose a “stencil”, based on r cells to the 
left, s cells to the right, and I t itself if r, s > 0, with r -f s + 1 = k: 

(2.8) ^(0 = {A — ri f»+a} • 

There is a unique polynomial of degree at most k - 1 = r + s, denoted by p(x) (we will drop the subscript 
i when it does not cause confusion), whose cell average in each of the cells in S(i ) agrees with that of v(x): 



This polynomial p(x) is the k - th order approximation we are looking for, as it is easy to prove (2.5), see the 
discussion below, as long sis the function v(x) is smooth in the region covered by the stencil S(i). 

For solving Problem 2.1, we also need the approximations to the values of v(x) at the cell boundaries, 

(2.6) . Since the mappings from the given cell averages vj in the stencil S(i) to the values vj and v U in 

(2.6) are linear, there exist constants cvj and Crj } which depend on the left shift of the stencil r of the stencil 
S(i) in (2.8), on the order of accuracy and on the cell sizes A xj in the stencil S*, but not on the function 
v itself, such that 

/c-l k - 1 

(2.10) ~ ^ , CrjVj— r+ji V i —2 = ^ 

3 j-o 2 j = 0 

We note that the difference between the values with superscripts ± at the same location a? i+ i is due to the 
possibility of different stencils for cell Ii and for cell /*+ j. If we identify the left shift r not with the cell I{ 
but with the point of reconstruction i.e. using the stencil (2.8) to approximate x i+ i, then we can drop 

the superscripts i and also eliminate the need to consider Crj in (2.10), as it is clear that 


Crj — Cr—ij. 

We summarize this as follows: given the k cell averages 


r» ***i — r+k — lj 

there are constants ev 3 such that the reconstructed value at the cell boundary 1 : 

fc-i 

(2.11) v i+h ~ y^Crj v i-r+j» 

i=0 
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is k - th order accurate: 


(2.12) v i+ ,=v(x i+h ) + 0(Ax k ). 

To understand how the constants {c^} are obtained, as well as how the accuracy property (2.5) is 
proven, we look at the primitive function of v(x): 

(2.13) V(x)= f v(Z)<%, 

J — OO 

where the lower limit — oo is not important and can be replaced by any fixed number. Clearly, V(x i+ i) can 
be expressed by the cell averages of v(x) using (2.4): 

(2.14) n* i+ i)= E r +iv (OdZ = E VjAxj, 

j--oo Jx j~^ j— OO 

thus with the knowledge of the cell averages {F;} we also know the primitive function V(x) at the cell 
boundaries exactly. If we denote the unique polynomial of degree at most A;, which interpolates V r (x^ + i) at 
the following k + 1 points: 


(2.15) 

by P(x), and denote its derivative by p(x): 


(2.16) 


p(x) = P'{x) , 


then it is easy to verify (2.9): 


1 

A xj 



* 2 

1 [ x i+i 

^7 J v{£)d£ = Vj, j = i-r,...,i + s, 


where the third equality holds because P(x) interpolates V(x) at the points and x J+ i whenever 

j = i — r, ...,i + s. This implies that p(x) is the polynomial we are looking for. Standard approximation 
theory (see an elementary numerical analysis book) tells us that 


P\x) = V'(x) + 0(Ax k ), X e Ii. 


This is the accuracy requirement (2.5). 

Now let us look at the practical issue of how to obtain the constants {crj} in (2.11). For this we could 
use the Lagrange form of the interpolation polynomial: 


(2.17) 


k 

^) = E^(wi) 

m~ 0 


n 

1 = 0 


x i- r+m-i ~ x i-r+l-$ 


l ^ m 
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For easier manipulation we subtract a constant V(x i _ r _ i) from (2.17), and use the fact that 

k k 


s n x 


x ~ x i-r+l-i 


m = 0 


1 = 0 
l ± m 


i— r+m-£ x i—r+l—% 


= 1 , 


to obtain: 
(2.18) 


P(x)-V(x i _ r _ i ) 


— X (^( x i-r+m-£) ^( x i-r-§)) IT x 

A 1 


x x i-r+l—\ 


m=0 


1 = 0 
l ^ m 


r+m— § x i— r+f— ^ 


Taking derivative on both sides of (2.18), and noticing that 


m — 1 


Vfo-r+m-j) - ^( x «-r-l) = X v i-r+j& x i-r+j 

j = 0 


because of (2.14), we obtain 


k m— 1 

(2.19) P{ x ) ” ^ ^ ^ ^ AXi_ r -f j 

m— 0 j=0 


(^ = 0 <7 = 0 (* x< - r+ «-0 


l ^ m q ^ mj 


n / = () ( X i-r+m-l X i-r+l -5) 

y l ^ m 


Evaluating the expression (2.19) at x = x iH _i , we finally obtain 


v i+ i =p(x i+ i) 


fc-l 

J=0 


( X *+§ x i-r+g-i) 


( yk prfc 

^ Z = 0 11 9 = 0 

m==j+l n I _ Q ( x i-r+m-J ~ x «-r+i— 


V 


Z = 0 
Z 7^ m 


^ x i - r 4 7 t — r f 7 i 


i.e. the constants in (2.11) are given by 


( = o 0 ( X i+i X i-r+,-i)) 


(2.20) 


Crj = 


Z = 0 -- q 

k l 7 ^ 77i q ^ m,l 




^2 

n=j+l n ^ _ Q £ — r+l— § ) 

l ^ m 


Ax*_ r _f-j . 


Although there are many zero terms in the inner sum of (2.20) when x i+ £ is a node in the interpolation, 
we will keep this general form so that it applies also to the case where x i+ i is not an interpolation point. 
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For a nonuniform grid, one would want to pre-compute the constants {c^} as in (2.20), for 0 <i<N, 
— I <r <k — 1, and 0 < j < k - 1, and store them before solving the PDE. 

For a uniform grid, Ax* = Ax, the expression for Crj does not depend on i or Ax any more: 


(2.21) 


°rj= £ 

m=j+l 


1 = 0 

1 1 ^ - m 


n q=0 (.<■-„+!) 


q ^ m,l 


if 


l = 0 

l ^ 771 


(m-l) 


We list in Table 2.1 the constants in this uniform grid case (2.21), for order of accuracy between 
k = 1 and k — 7. 

From Table 2.1, we would know, for example, that 


1„ 5_ 1_ 

+ 3 V *+i + 


0( Ax 3 ) . 


2.1.2. Conservative approximation to the derivative from point values. The second approx- 
imation problem we will face, in solving hyperbolic conservation laws using point values (finite difference 
schemes, see Sect. 2.3.2), is the following problem in obtaining high order conservative approximation to the 
derivative from point values [69, 70]. 


Problem 2.2. One dimensional conservative approximation. 

Given the point values of a function v(x): 

(2.22) Vi = v(xi), i = 1,2,..., AT, 

find a numerical flux function 


(2.23) 


v i+i 


= v(Vi- r 


u i + s 


), i = 0,l,...,N , 


such that the flux difference approximates the derivative v l (x) to k - th order accuracy: 

(2.24) (v <+ i -Vi) =v'(x0 + O(Ax fc ), i = 0, 1, 

We again ignore the boundary conditions here and assume that Vi is available for i < 0 and i > N if 
needed. 

The solution of this problem is essential for the high order conservative schemes based on point values 
(finite difference) rather than on cell averages (finite volume). 

This problem looks quite different from Problem 2.1. However, we will see that there is a close relationship 
between these two. We assume that the grid is uniform , Ax* = Ax. This assumption is, unfortunately, 
essential in the following development. 

If we can find a function h(x), which may depend on the grid size Ax, such that 


(2.25) 
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then clearly 



hence all we need to do is to use 

(2.26) 6 i+i =M*, + i) + 0(A* fc ) 

to achieve (2.24). We note here that it would look like an 0( Ax* +1 ) term in (2.26) is needed in order to get 
(2.24), due to the Ax term in the denominator. However, in practice, the 0(Ax*) term in (2.26) is usually 
smooth, hence the difference in (2.24) would give an extra 0(Ax), just to cancel the one in the denominator. 

It is not easy to approximate h(x) via (2.25), as it is only implicitly defined there. However, we notice 
that the known function v(x) is the cell average of the unknown function h(x ), so to find h(x) we just need 
to use the reconstruction procedure described in Sect. 2.1.1. If we take the primitive of h{x)\ 

(2.27) H(x)= f X h(Odl £, 

J — oc 

then (2.25) clearly implies 

(2.28) H(x i+i ) = J2 MfR = Ax £ v 3 . 

j~-OC x j - £ j — oo 

Thus, given the point values we “identify” them as cell averages of another function h(x) in (2.25), 
then the primitive function H (x) is exactly known at the cell interfaces x = x i+ 1 . We thus use the same 
reconstruction procedure described in Sect. 2.1,1, to get a fc-th order approximation to /i(x 1+ i), which is 
then taken as the numerical flux t) i+ i in (2.23). 

In other words, if the “stencil” for the flux in (2.23) is the following k points: 

(2.29) Xt_ r , 3?i+5 > 
where r + s = k — 1, then the flux Vi+i is expressed as 

k-l 

( 2 - 30 ) 5 i+ %=Yl l Cr i Vi - r +j' 

j - 0 

where the constants {crj} are given by (2.21) and Table 2.1. 

From Table 2.1 we would know, for example, that if 

1 5 1 

u i-f i = + 3^+1 , 

then 

^ (t) j+ ! - ®i_j) = v'(xi) + 0( Ax 3 ). 

We emphasize again that, unlike in the reconstruction procedure in Sect. 2.1.1, here the grid must be 
uniform: A Xj = Ax. Otherwise, it can be proven that no choice of constants in (2.30) (which may 
depend on the local grid sizes but not on the function v(x)) could make the conservative approximation to 
the derivative (2.24) higher than second order accurate ( k > 2). The proof is a simple exercise of Taylor 
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expansions. Thus, the high order finite difference (third order and higher) discussed in these lecture notes 
can apply only to uniform or smoothly varying grids. 

Because of this equivalence of obtaining a conservative approximation to the derivative (2. 23)- (2. 24) and 
the reconstruction problem discussed in Sect. 2.1.1, we will only need to consider the reconstruction problem 
in the following sections. 


2.1.3. Fixed stencil approximation. By fixed stencil, we mean that the left shift r in (2.8) or (2.29) 
is the same for all locations i. Usually, for a globally smooth function v(x), the best approximation is 
obtained either by a central approximation r = s — 1 for even k (here central is relative to the location x i+ 1 ), 
or by a one point upwind biased approximation r = s or r = s — 2 for odd k. For example, if the grid is 
uniform Axi = Ax, then a central 4th order reconstruction for t^+i, in (2.11), is given by 


1 _ 7 _ 7 _ 1 _ ^ /A 

i + \2 Vi + 12 ~ \2 Vi+2 + °^ AX ^ ’ 

and the two one point upwind biased 3rd order reconstructions for in (2.11), are given by 

1- 5_ 1_ ^ /A 

- i + -Vi + “U i+ i + 0(Ax ) 

1_ 5_ 1_ ^ /A 

or v i+ i= -Vi + -v i+1 - -v i+2 + 0{ Ax ) . 


Similarly, a central 4th order flux (2.30) is 


= -^-1 + —Vi + — v i+1 


12 


12 


12 


12 


Vt+2 : 


which gives 

^ - Vi) = v '( x i) + 0{ Ax 4 ), 

and the two one point upwind biased 3rd order fluxes (2.30) are given by 


or 


1 5 1 

v *+ J = + Q v i + 3 V i+i 

1 5 1 


which gives 

2^ ( € *+^ “ °i-i) = v '( x i) + 0(Ax 3 ). 

Traditional central and upwind schemes, either finite volume or finite difference, can be derived by these 
fixed stencil reconstructions or flux differenced approximations to the derivatives. 


2.2. ENO and WENO Approximations in ID. For solving hyperbolic conservation laws, we are 
interested in the class of piecewise smooth functions. A piecewise smooth function v(x) is smooth (i.e. it 
has as many derivatives as the scheme calls for) except for at finitely many isolated points. At these points, 
v(x) and its derivatives are assumed to have finite left and right limits. Such functions are “generic” for 
solutions to hyperbolic conservation laws. 

For such piecewise smooth functions, the order of accuracy we refer to in these lecture notes are formal , 
that is, it is defined as whatever accuracy determined by the local truncation error in the smooth regions of 
the function. 
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Fig. 2.1. Fixed central stencil cubic interpolation (left) and ENO cubic interpolation (right) for the step function. Solid: 
exact function; Dashed: interpolant piecewise cubic polynomials. 


If the function v(x) is only piecewise smooth, a fixed stencil approximation described in Sect. 2.1.3 
may not be adequate near discontinuities. Fig. 2.1 (left) gives the 4-th order (piecewise cubic) interpolation 
with a central stencil for the step function, i.e. the polynomial approximation inside the interval [ar^i, x i+ i] 
interpolates the step function at the four points | , x i _ i , x i+ 1 , x i+ 3 . Notice the obvious over /undershoots 

for the cells near the discontinuity. 

These oscillations (termed the Gibbs Phenomena in spectral methods) happen because the stencils, as 
defined by (2.15), actually contain the discontinuous cell for Xi close enough to the discontinuity. As a result, 
the approximation property (2.5) is no longer valid in such stencils. 

2.2.1. ENO approximation. A closer look at Fig. 2.1 (left) motivates the idea of “adaptive stencil”, 
namely, the left shift r changes with the location Xi. The basic idea is to avoid including the discontinuous 
cell in the stencil, if possible. 

To achieve this effect, we need to look at the Newton formulation of the interpolation polynomial. 

We first review the definition of the Newton divided differences. The 0-th degree divided differences of 
the function V(x) in (2.13)-(2.14) are defined by: 


(2-31) 

and in general the j-th degree divided differences, for j > 1, are defined inductively by 

(2.32) V\x,_ h ..^,,1- 

Similarly, the divided differences of the cell averages v in (2.4) are defined by 
(2-33) *[*<]=% 


and in general 


(2.34) 




v\Xi ± 1 


We note that, by (2.14), 


Xj+j\ v\Xj) ..., 
X{+j X{ 


(2.35) 


»*«+*] = 


V(x i+ r) -Vfoi- U 

X i+ \ ~ X i-£ 


= Vi, 
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i.e. the O-th degree divided differences of v are the first degree divided differences of V{x). We can thus 
write the divided differences of V(x) of first degree and higher by those of v of 0-th degree and higher, using 
(2.35) and (2.32). 

The Newton form of the fc-th degree interpolation polynomial P(x), which interpolates V(x) at the k- f 1 
points (2.15), can be expressed using the divided differences (2.31)-(2.32) by 


(2.36) 


k 


P(x) = 

j=0 


x i-r+j- 


ii n (* 

m=0 





We can take the derivative of (2.36) to get p(x) in (2.16): 


(2.37) 


p( x ) = 


j=l 


J- 1 

J-l 

n 

m=0 

O 

II 



i—r+l— 



l / m 


Notice that only first and higher degree divided differences of V{x) appear in (2.37). Hence by (2.35), we 
can express p(x) completely by the divided differences of v, without any need to reference V (x) . 

Let us now recall an important property of divided differences: 

(2.38) V\xi_ 7 

for some £ inside the stencil: x,_i < £ < x t ' +J -_ i , as long as the function V(x) is smooth in this stencil. If 
V(x) is discontinuous at some point inside the stencil, then it is easy to verify that 

(2.39) y [xi-i, - >Zi + ;-}] = C> (^j) • 

Thus the divided difference is a measurement of the smoothness of the function inside the stencil. 

We now describe the ENO idea by using (2.36). Suppose our job is to find a stencil of k + 1 consecutive 
points, which must include x t _ i and x i+ j, such that V(x) is “the smoothest” in this stencil comparing with 
other possible stencils. We perform this job by breaking it into steps, in each step we only add one point to 
the stencil. We thus start with the two point stencil 


(2.40) ^ 2 (*) — { x i- j ) X i+% }> 

where we have used S to denote a stencil for the primitive function V . Notice that the stencil S for V has 
a corresponding stencil S for v through (2.35), for example (2.40) corresponds to a single cell stencil 


5(0 = {/.} 

for v. The linear interpolation on the stencil S 2 (i) in (2.40) can be written in the Newton form as 

P\x) = V[xi_£,x j+i ] (x-x^i) . 

At the next step, we have only two choices to expand the stencil by adding one point: we can either add the 
left neighbor resulting in the following quadratic interpolation 

(2.41) R(x) = P 1 (x) + V [x 4 _ 3 , Xi_ 1 , x i+ 1 ] (x - x t _ 1 ) (x - x j+ i ) , 
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or add the right neighbor x i+ s , resulting in the following quadratic interpolation 
(2.42) S(x) = P (x) + V[x ^_^ , 1 , x i+ |] ^x — — x i+ i^ * 

We note that the deviations from P l (x) in (2.41) and (2.42), are the same function 



multiplied by two different constants 


(2.43) 


1 5 X % — i J X 




and V\xi_ i , > ^*+ §]* 


These two constants are the two second degree divided differences of V'(x) in two different stencils. We 
have already noticed before, in (2.38) and (2.39), that a smaller divided difference implies the function is 
“smoother” in that stencil. We thus decide upon which point to add to the stencil, by comparing the two 
relevant divided differences (2.43), and picking the one with a smaller absolute value. Thus, if 


(2.44) 


n 


Xi_ 




we will take the 3 point stencil as 


^ 3(0 — 

otherwise, we will take 


^3(0 — i X i-% > J }* 

This procedure can be continued, with one point added to the stencil at each step, according to the 
smaller of the absolute values of the two relevant divided differences, until the desired number of points in 
the stencil is reached. 

We note that, for the uniform grid case Ax* = Ax, there is no need to compute the divided differences 
as in (2.32). We should use undivided differences instead: 

(2.45) V >=V[x i _i,x i+ i] = v i 

(see (2.35)), and 

V < 1 , ..., Zj+j+i > = V < ^t+j+i > -V < >, 

(2.46) j > 1 . 

The Newton interpolation formulae (2.36)- (2.37) should also be adjusted accordingly. This both saves com- 
putational time and reduces round-off effects. 

The FORTRAN program for this ENO choosing process is very simple: 

* assuming the m-th degree divided (or undivided) differences 

* of V(x), with x_i as the left-most point in the arguments, 

* axe stored in V(i,m), also assuming that "is" is the 

* left-most point in the stencil for cell i for a k-th 

* degree polynomial 
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is=i 

do m=2,k 

if (abs (V (is-1 ,m) ) . It . abs (V(is ,m) ) ) is=is-l 

enddo 

Once the stencil S(i) } hence 5(i), in (2.8) is found, one could use (2.11), with the prestored values of 
the constants , (2.20) or (2.21), to compute the reconstructed values at the cell boundary. Or, one could 
use (2.30) to compute the fluxes. An alternative way is to compute the values or fluxes using the Newton 
form (2.37) directly. The computational cost is about the same. 

We summarize the ENO reconstruction procedure in the following 

Procedure 2.1. ID ENO reconstruction* 

Given the cell averages {v*} of a function u(x), we obtain a piecewise polynomial reconstruction, of 
degree at most k - 1, using ENO, in the following way: 

1. Compute the divided differences of the primitive function V(x), for degrees 1 to k , using v , (2.35) 
and (2.32). 

If the grid is uniform Ax t = Ax, at this stage, undivided differences (2.45)-(2.46) should be computed 
instead. 

2. In cell I i} start with a two point stencil 

5 2 (i) = 

for V'(x), which is equivalent to a one point stencil, 


for v . 

3. For l = 2, ..., k , assuming 


Si(i) — {Xj+i > Xj + j_i } 


is known, add one of the two neighboring points, x^i or x j+l +i } to the stencil, following the ENO 
procedure: 

• If 


(2.47) 




add to the stencil Si(i) to obtain 

s+i(0 = j x j+i— ^ }> 


• Otherwise, add x J+ j + i to the stencil Si(i) to obtain 

Si+i{i) = {x j+ i,...,x J+ | + i}. 


4. Use the Lagrange form (2.19) or the Newton form (2.37) to obtain pi(x), which is a polynomial of 
degree at most k - 1 in /*, satisfying the accuracy condition (2.5), as long as v(x) is smooth in 
We could use Pi{x) to get the approximations at the cell boundaries: 


« i+ l = PiK+l). v t-i = 
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However, it is usually more convenient, when the stencil is known, to use (2.10), with Cr 3 defined 
by (2.20) for a nonuniform grid, or by (2.21) and Table 2.1 for a uniform grid, to compute an 
approximation to v(x ) at the cell boundaries. 

For the same piecewise cubic interpolation to the step function, but this time using the ENO procedure 
with a two point stencil § 2 (i) = in the Step 2 of Procedure 2.1, we obtain a non-oscillatory 

interpolation, in Fig. 2.1 (right). 

For a piecewise smooth function F(x), ENO interpolation starting with a two point stencil S 2 (i) = 
{Zi_i , £*+ 1 } in the Step 2 of Procedure 2.1, as was shown in Fig. 2.1 (right), has the following properties 
[39]: 2 

1. The accuracy condition 

Pi{x) = V{x) + 0(Ax fc+1 ), X € Ii 
is valid for any cell Ii which does not contain a discontinuity. 

This implies that the ENO interpolation procedure can recover the full high order accuracy right up 
to the discontinuity. 

2. Pi(x) is monotone in any cell Ii which does contain a discontinuity of V(x). 

3. The reconstruction is TVB (total variation bounded). That is, there exists a function z(x), satisfying 

z(x) = Pi(x) + 0( Ax k+1 ), x € Ii 

for any cell /{, including those cells which contain discontinuities, such that 


TV(z) <TV(V). 


Property 3 is clearly a consequence of Properties 1 and 2 (just take z(x) to be V(x) in the smooth cells 
and take z{x) to be Pi(x) in the cells containing discontinuities). It is quite interesting that Property 2 holds. 
One would have expected trouble in those “shocked cells”, i.e. cells Ii which contain discontinuities, for ENO 
would not help for such cases as the stencil starts with two points already containing a discontinuity. We 
will give a proof of Property 2 for a simple but illustrative case, i.e. when V (:r) is a step function 


V(x) = 


{ 


0 , 

1 , 


x < 0; 
x > 0. 


and the k-th degree polynomial P(x) interpolates V(x) at k + 1 points 


containing the discontinuity 


X 1 < x 3 < ... < XL.il 
2 3 2 


X io-i < 0 < X Jo+i 


for some jo between 1 and k. For any interval which does not contain the discontinuity 0: 


( 2 -48) j ^ j 0 , 

we have 


P(x M ) = V(x w ) = V(x J+i ) = P(x iH ), 
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hence there is at least one point in between, Xj_i < < x J+ i, such that P' (£j ) = 0. This way we can 

find k - 1 distinct zeroes for P'(x), as there are k - 1 intervals (2.48) which do not contain the discontinuity 

0. However, P'(x) is a non-zero polynomial of degree at most k - 1, hence can have at most k - 1 distinct 
zeroes. This implies that P'(x) does not have any zero inside the shocked interval [x jo _^ he. P(x) 

is monotone in this shocked interval. This proof can be generalized to a proof for Property 2 [39]. 

2.2.2. WENO approximation. In this subsection we describe the recently developed WENO (weighted 
ENO) reconstruction procedure [53, 43]. WENO is based on ENO, of course. For simplicity of presentation, 
in this subsection we assume the grid is uniform, i.e. Ax{ = Ax. 

As we can see from Sect. 2.2.1, ENO reconstruction is uniformly high order accurate right up to the 
discontinuity. It achieves this effect by adaptively choosing the stencil based on the absolute values of divided 
differences. However, one could make the following remarks about ENO reconstruction, indicating rooms 
for improvements: 

1. The stencil might change even by a round-off error perturbation near zeroes of the solution and its 
derivatives. That is, when both sides of (2.47) are near 0, a small change at the round off level 
would change the direction of the inequality and hence the stencil. In smooth regions, this “free 
adaptation” of stencils is clearly not necessary. Moreover, this may cause loss of accuracy when 
applied to a hyperbolic PDE [63, 67]. 

2. The resulting numerical flux (2.23) is not smooth, as the stencil pattern may change at neighboring 
points. 

3. In the stencil choosing process, k candidate stencils are considered, covering 2k - 1 cells, but only 
one of the stencils is actually used in forming the reconstruction (2.10) or the flux (2.30), resulting in 
Jfc-th order accuracy. If all the 2k - 1 cells in the potential stencils are used, one could get (2k - l)-th 
order accuracy in smooth regions. 

4. ENO stencil choosing procedure involves many logical “if’ structures, or equivalent mathematical 
formulae, which are not very efficient on certain vector computers such as CRAYs (however they are 
friendly to parallel computers). 

There have been attempts in the literature to remedy the first problem, the “free adaptation” of stencils. 
In [28] and [67], the following “biasing” strategy was proposed. One first identity a “preferred” stencil 


(2.49) 


Spref{i) — ...» X^_ r _|_ fe+ 1 } , 


which might be central or one-point upwind. One then replaces (2.47) by 


V\x, 


j- £>*••> X j+l- 


< b |f[x j+ i, ...,x^ +i+ i] , 


if 


x j+ \ > X i-r+% > 

i.e. if the left- most point in the current stencil Si(i) has not reached the left- most point x i _ r+ x of the 

preferred stencil 5 , pT - c /(2) in (2.49) yet; otherwise, if 


— x »-r+i 3 


one replaces (2.47) by 


aw ”* x j+i - •••> | • 
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Here, b > 1 is the so-called biasing parameter. Analysis in [67] indicates a good choice of the parameter 
b = 2. The philosophy is to stay as close as possible to the preferred stencil, unless the alternative candidate 
is, roughly speaking, a factor b > 1 better in smoothness. 

WENO is a more recent attempt to improve upon ENO in these four points. The basic idea is the 
following: instead of using only one of the candidate stencils to form the reconstruction, one uses a convex 
combination of all of them. To be more precise, suppose the k candidate stencils 

( 2 - 50 ) S r (i) = {Xi-r^.^Xi-r+k-i}, r = 0 1 


produce k different reconstructions to the value u i+ ±, according to (2.11), 


k - 1 


(2.51) 


•+A = E CrjVi-r+j, T = 0 , k - 1 , 


j= 0 


WENO reconstruction would take a convex combination of all defined in (2.51) as a new approximation 

1 ' 3 

to the cell boundary value u(x i+ ±): 




(2.52) 




r = 0 


Apparently, the key to the success of WENO would be the choice of the weights uj t . We require 


k - 1 


(2.53) 


uj r ^ 0 , ujf — 1 


r=0 


for stability and consistency. 

If the function v(x) is smooth in all of the candidate stencils (2.50), there arc constants d r such that 

(2.54) v iH = = v(x l+h )+0(Ax 2k -'). 

r=0 

For example, d T for 1 < k < 3 are given by 

do = 1 , k z= 1 j 
J 2 J 1 
do 2 ’ ™ — 2 i 

do= 10’ dl = 5’ d2 = To’ fc = 3 

We can see that d r is always positive and, due to consistency, 

k - 1 


(2.55) 


d r — i . 


r = 0 


In this smooth case, we would hke to have 
(2.56) a > r = d r + 0{ Ax*' 1 ), r = 0, ..., k - 1 , 


which would imply (2k — l)-th order accuracy: 


(2.57) 


k-i 


V i+ x = 


r= 0 


v(x i+ i) + 0(Ax 2k *) 
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because 


k - 1 


E 


k - 1 


UrV l+i 


W. 


r=0 


i+J 


fc-1 


= ^(w r 



A: — 1 

= ^ 0(Ax fe_1 ) 0(Ax k ) = 0(Ax 2fc_1 ) 

r = 0 


where in the first equality we used (2.53) and (2.55). 

When the function v(x) has a discontinuity in one or more of the stencils (2.50), we would hope the 
corresponding weight (s) u r to be essentially 0, to emulate the successful ENO idea. 

Another consideration is that the weights should be smooth functions of the cell averages involved. In 
fact, the weights designed in [43] and described below are C°° . 

Finally, we would like to have weights which are computationally efficient. Thus, polynomials or rational 
functions are preferred over exponential type functions. 

All these considerations lead to the following form of weights: 


(2.58) 


UJ r = 


a r 


E k —1 
3=0 


a s 


r = 0 , k — 1 


with 

(2.59) 


dr 

“ r= («:+&)*■ 


Here e > 0 is introduced to avoid the denominator to become 0. We take e = 10 6 in all our numerical 
tests [43]. p r are the so-called “smooth indicators” of the stencil S r (i): if the function v(x) is smooth in the 
stencil S r (t), then 


p r = 0( Ax 2 ), 


but if r»(x) has a discontinuity inside the stencil S r ( t), then 

Pr = 0(1). 


Translating into the weights u T in (2.58), we will have 

u) r — 0 ( 1 ) 

when the function u(x) is smooth in the stencil S r (i), and 

U r — 0( Ax 4 ) 

if v(x) has a discontinuity inside the stencil 5 r (z). Emulation of ENO near a discontinuity is thus achieved. 

One also has to worry about the accuracy requirement (2.56), which must be checked when the specific 
form of the smooth indicator /3 r is given. For any smooth indicator /3 r , it is easy to see that the weights 
defined by (2.58) satisfies (2.53). To satisfy (2.56), it suffices to have, through a Taylor expansion analysis: 

(2.60) Pr = D( 1 + 0( Ax*' 1 )), r = 0, ..., k - 1, 

where D is a nonzero quantity independent of r (but may depend on Ax). 
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As we have seen in Sect. 2.2.1, the ENO reconstruction procedure chooses the “smoothest” stencil by 
comparing a hierarchy of divided or undivided differences. This is because these differences can be used 
to measure the smoothness of the function on a stencil, (2. 38)- (2. 39). In [43], after extensive experiments, 
a robust (for third and fifth order at least) choice of smooth indicators (3 r is given. As we know, on each 
stencil S r (i), we can construct a (k — l)-th degree reconstruction polynomial, which if evaluated at x = x i+ i, 
renders the approximation to the value v(x i+ 1 ) in (2.51). Since the total variation is a good measurement for 
smoothness, it would be desirable to minimize the total variation for this reconstruction polynomial inside 
Ii. Consideration for a smooth flux and for the role of higher order variations leads us to the following 
measurement for smoothness: let the reconstruction polynomial on the stencil S r (i) be denoted by p r (x), we 
define 


*-i rx 

-E/ 

[-1 Jx i- 


Ax 21 - 1 ( d* Pr(x) 
-* V d‘x 


The right hand side of (2.61) is just a sum of the squares of scaled L 2 norms for all the derivatives of the 
interpolation polynomial p r (x) over the interval (x^i, x i+ i). The factor Ax 2 *" 1 is introduced to remove 
any Ax dependency in the derivatives, in order to preserve self-similarity when used to hyperbolic PDEs 
(Sect. 2.3). 

We remark that (2.61) is similar to but smoother than the total variation measurement based on the L 1 
norm. It also renders a more accurate WENO scheme for the case k = 2 and 3. 

When k — 2, (2.61) gives the following smoothness measurement [53, 43]: 

0o = (Vi+1 - Vi) 2 , 

(2.62) /?i = (t'i — Vi-i) 2 . 

For k = 3, (2.61) gives [43]: 

13 1 

0° = - 2rJ i+ i + v i+2 ) 2 + -(3Fj - 4v i+ i + v i+2 ) 2 , 

( 2 - 63 ) Pi = ~ 2tJj + Vi + i) 2 + -Ui+i) 2 , 

13 1 

02 = — (Vi — 2 - 2Ui_i + Vi) 2 + -{Vi - 2 - 4Vj_i + 3t7j) 2 . 

We can easily verify that the accuracy condition (2.60) is satisfied, even near smooth extrema [43]. This 
indicates that (2.62) gives a third order WENO scheme, and (2.63) gives a fifth order one. 

Notice that the discussion here has a one point upwind bias in the optimal linear stencil, suitable for 
a problem with wind blowing from left to right. If the wind blows the other way, the procedure should be 
modified symmetrically with respect to x i+ i. 

In summary, we have the following WENO reconstruction procedure: 

Procedure 2.2. ID WENO reconstruction. 

Given the cell averages {Fj} of a function v(x), for each cell 7^, we obtain upwind biased (2k — l)-th 
order approximations to the function v(x) at the cell boundaries, denoted by A and v7 2 , in the following 
way: 

1. Obtain the k reconstructed values of &-th order accuracy, in (2.51), based on the stencils 

(2.50), for r — 0, ..., k - 1; 

Also obtain the k reconstructed values of A;-th order accuracy, using (2.10), again based on 

the stencils (2.50), for r = 0, ..., k — 1; 
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2. Find the constants d r and <f r , such that (2.54) and 

Vi-$ = Y; dr v \ r J . i = + 0(Ax 2fc_1 ) 

r= 0 3 

are valid. By symmetry, 


d r — dk—i— r . 


3. Find the smooth indicators /? r in (2.61), for all r = 0 - 1. Explicit formulae for A: = 2 and 
A; = 3 are given in (2.62) and (2.63) respectively. 

4. Form the weights uj t and u T using (2. 58)- (2. 59) and 


a r 


U) r = 


2 ^ 5=0 


a r = 


d T 


(e + /3 r ) 2 ' 


r = 0, k — 


5. Find the (2k — l)-th order reconstruction 

fc-i fc-i 

(2.64) v~ + ^ u ' v i+k > v t J = E 

r=0 


r— 0 


We can obtain weights for higher orders of k (corresponding to seventh and higher order WENO schemes) 
using the same recipe. However, these schemes of seventh and higher order have not been extensively tested 
yet. 

2.3. ENO and WENO Schemes for ID Conservation Laws. In this section we describe the ENO 
and WENO schemes for ID conservation laws: 


(2.65) 


U t (x y t) + fx(u(x,t)) = 0 


equipped with suitable initial and boundary conditions. 

We will concentrate on the discussion of spatial discretization, and will leave the time variable t contin- 
uous (the method-of-lines approach). Time discretizations will be discussed in Sect. 4.2. 

Our computational domain is a < x < b. We have a grid defined by (2.1), with the notations (2. 2)- (2. 3). 
Except for in Sect. 2.3.3, we do not consider boundary conditions. We thus assume that the values of the 
numerical solution are also available outside the computational domain whenever they are needed. This 
would be the case for periodic or compactly supported problems. 

2.3.1. Finite volume formulation in the scalar case. For finite volume schemes, or schemes based 
on cell averages, we do not solve (2.65) directly, but its integrated version. We integrate (2.65) over the 
interval U to obtain 

(2-66) (/(«(x i+ i,t)-/(u(x f _ f ,t))) , 


where 

(2.67) 



is the cell average. We approximate (2.66) by the following conservative scheme 


( 2 . 68 ) 


du{(t ) 

dt 


1 

A Xi 


(/ i+ *-/i-i). 
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where Ui(t) is the numerical approximation to the cell average u(x,,t), and the numerical flux f i+ 1 is defined 
by 

(2 ’ 69) fi+l = h { u 7 + l’< + i) 

with the values uf + ^ obtained by the ENO reconstruction Procedure 2.1, or by the WENO reconstruction 
Procedure 2.2. 

The two argument function h in (2.69) is a monotone flux. It satisfies: 

• h(a , b ) is a Lipschitz continuous function in both arguments; 

• h{a,b) is a nondecreasing function in a and a nonincreasing function in b . Symbolically /i(f, j); 

• h(a,b) is consistent with the physical flux /, that is, h(a,a ) = /(a). 

Examples of monotone fluxes include: 

1. Godunov flux: 


(2.70) 


2. Engquist-Osher flux: 


h(a, 6) 


mm a < u < b f(u ) if a < b 

maxfc< u < a f(u) if a > b 


(2.71) 


h(a, 6) = / max(/'(u), 0)du + 

J o 



min(/'(u),0)dif + /( 0). 


3. Lax- Friedrichs flux: 


(2.72) h(a, b) = i [/(a) + /(&) - a(6 - a)] 

where a = max u |/'(u)| is a constant. The maximum is taken over the relevant range of u. 

We have listed the monotone fluxes from the least dissipative (less smearing of discontinuities) to the most. 
For lower order methods (order of reconstruction is 1 or 2), there is a big difference between results obtained 
by different monotone fluxes. However, this difference becomes much smaller for higher order reconstructions. 
In Fig. 2.2, we plot the results of a right moving shock for the Burgers 1 equation ( f(u ) = in (2.65)), 
with first order reconstruction using Godunov and Lax-Friedrichs monotone fluxes (top), and with fourth 
order ENO reconstruction using Godunov and Lax- Friedrichs monotone fluxes (bottom). We can clearly sec 
that, while the Godunov flux behaves much better for the first order scheme, the two fourth order ENO 
schemes behave similarly. We thus use the simple and inexpensive Lax-Friedrichs flux in most of our high 
order calculations. 

We remark that, by the classic Lax-Wendroff theorem [51], the solution to the conservative scheme 
(2.68), if converges , will converge to a weak solution of (2.65). 

In summary, to build a finite volume ENO scheme (2.68), given the cell averages {fZ*} (we will often 
drop the explicit reference to the time variable t) y we proceed as follows: 


Procedure 2.3. Finite volume ID scalar ENO and WENO. 

1. Follow the Procedure 2.1 in Sect. 2.2.1 for ENO, or the Procedure 2.2 in Sect. 2.2.2 for WENO, to 

obtain the k - th order reconstructed values u7 A and j for all i ; 

a *+5 

2. Choose a monotone flux (e.g., one of (2.70) to (2.72)), and use (2.69) to compute the flux for 
all i; 

3. Form the scheme (2.68). 

Notice that the finite volume scheme can be applied to arbitrary nonuniform grids. 
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Fig. 2.2. First order (top) and fourth order (bottom) ENO schemes for the Burgers equation, with the Godunov flux (left) 
and the Lax- Friedrichs flux (right). Solid lines: exact solution; Circles: the computed solution at t = 4. 


2.3.2. Finite difference formulation in the scalar case. We first assume the grid is uniform and 
solve (2.65) directly using a conservative approximation to the spatial derivative: 

dui(t) 1 (} i \ 

(2.73) — = 

where Ui(t) is the numerical approximation to the point value u{xi,t), and the numerical flux 

/i-f 1 = f { u i-ri 


satisfies the following conditions: 

• / is a Lipschitz continuous function in all the arguments; 

• f is consistent with the physical flux /, that is, /(u, = f(u). 

Again the Lax-Wendroff theorem [51] applies. The solution to the conservative scheme (2.73), if con- 
verges , will converge to a weak solution of (2.65). 

The numerical flux / i+A is obtained by the ENO or WENO reconstruction procedures, Procedure 2.1 
or 2.2, with v{x) = /(u(x,t)). For stability, it is important that upwinding is used in constructing the flux. 
The easiest and the most inexpensive way to achieve upwinding is the following: compute the Roe speed 


(2.74) 


a i+i = 


f{ u i+ 1) f{ u i) 

tti+i - Ui 


• if a t+ ^ > 0, then the the wind blows from the left to the right. We would use v i j r i f° r the numerical 
flux / i+i ; 
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• if a i+i < 0, then the wind blows from the right to the left. We would use 
flux f i+ 1 . 


for the numerical 


This produces the Roe scheme [62] at the first order level. For this reason, the ENO scheme based on this 
approach was termed “ENO-Roe” in [70]. 

In summary, to build a finite difference ENO scheme (2.73) using the ENO-Roe approach , given the point 
values {m} (we again drop the explicit reference to the time variable £), we proceed as follows: 


Procedure 2.4. Finite difference ID scalar ENO- and WENO-Roe. 

1. Compute the Roe speed a i+ i for all i using (2.74); 

2. Identify Vi = f(ui ) and use the ENO reconstruction Procedure 2.1 or the WENO reconstruction 
Procedure 2.2, to obtain the cell boundary values if a i+ i > 0, or if < 0; 

3. If the Roe speed at x i+ 1 is positive 

a i+ i > 0, 

then take the numerical flux as: 




otherwise, take the the numerical flux as: 


fi + i = <+i' 

4. Form the scheme (2.73). 

One disadvantage of the ENO-Roe approach is that entropy violating solutions may be obtained, just 
like in the first order Roe scheme case. For example, if ENO-Roe is applied to the Burgers equation 


u t + 


(*).- 


with the following initial condition 


u(x , 0) 


— 1, if x < 0, 
1, if x > 0, 


it will converge to the entropy violating expansion shock: 


u(x , t) 


— 1, if x < 0, 
1, if x > 0. 


Local entropy correction could be used to remedy this [70]. However, it is usually more robust to use a 
global “flux splitting” : 


(2.75) 


f(u) = /+(«) + / (u) 


where 


(2.76) 


rf / + ( M ) > n d f ( M ) ^ n 

du ~ ’ du ~ 


We would need the positive and negative fluxes f ± ( u ) to have as many derivatives as the order of the scheme. 
This unfortunately rules out many popular fluxes (such as those of van Leer [79] and Osher [58]) for high 
order methods in this framework. 
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The simplest smooth splitting is the Lax-Friedrichs splitting: 

(2.77) / ± (u) = ^(/( u ) ± au ) 

where a is again taken as a — max u |/'(u)| over the relevant range of u . 

We note that there is a close relationship between a flux splitting (2.75) and a monotone flux (2.69). In 
fact, for any flux splitting (2.75) satisfying (2.76), 

(2.78) h(ct,b) = / + (a) + /~(&) 

is clearly a monotone flux. However, not every monotone flux can be written in the flux split form (2.75). 
For example, the Godunov flux (2.70) cannot. 

With the flux splitting (2.75), we apply the the ENO or WENO reconstruction procedures, Procedure 
2.1 or 2.2, with v(x) = / + (u(x,t)) and u(z) = f~{u(x } t)) separately, to obtain two numerical fluxes ^ 

and ! , and then sum them to get the numerical flux / i+ 1 . 

In summary, to build a finite difference (FD) ENO or WENO scheme (2.73) using the flux splitting 
approach , given the point values we proceed as follows: 

Procedure 2.5. FD ID scalar flux splitting ENO and WENO. 

1. Find a smooth flux splitting (2.75), satisfying (2.76); 

2. Identify Vi = / + (u*) and use the ENO or WENO reconstruction procedure, Procedure 2.1 or 2.2, to 
obtain the cell boundary values for all z; 

3. Take the positive numerical flux as 

ft* = 

4. Identify l h = f~(u t ) and use the ENO or WENO reconstruction procedures, Procedure 2.1 or 2.2, 

to obtain the cell boundary values v + ' , for all i\ 

2 

5. Take the negative numerical flux as 

6. Form the numerical flux as 

= + h +^ 5 

7. Form the scheme (2.73). 

We remark that the finite difference scheme in this section and the finite volume scheme in Sect. 2.3.1 
are equivalent for one dimensional, linear PDE with constant coefficients: the only difference is in the initial 
conditions (one uses point values and the other uses cell averages of the exact initial condition). Notice 
that the schemes are still nonlinear in this case. However, this equivalence does not hold for a nonlinear 
PDE. Moreover, we will see later that there are significant differences in efficiency of the two approaches for 
multidimensional problems. 

In the following we test the accuracy of the fifth order finite difference WENO schemes on the linear 
equation: 

u t + u x = 0, — 1 < x < 1 

u(x, 0) = uq(x) periodic. 
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Table 2.2 

Accuracy on ut + u x = 0 with txo(^) = sin( 7 rar). 


Method 

N 

Loo error 

Loo order 

L\ error 

L\ order 


10 

2.98e-2 

- 

1.60e-2 

- 


20 

1.45e-3 

4.36 

7.41e-4 

4.43 

WENO- 5 

40 

4.58e-5 

4.99 

2.22e-5 

5.06 


80 

1.48e-6 

4.95 

6.91e-7 

5.01 


160 

4.41e-8 

5.07 

2.17e-8 

4.99 


320 

1.35e-9 

5.03 

6.79e-10 

5.00 


10 

4.98e-3 

- 

3.07e-3 

- 


20 

1.60e-4 

4.96 

9.92e-5 

4.95 

CENTRAL-5 

40 

5.03e-6 

4.99 

3.14e-6 

4.98 


80 

1.57e-7 

5.00 

9.90e-8 

4.99 


160 

4.91e~9 

5.00 

3.11e-9 

4.99 


320 

1.53e-10 

5.00 

9.73e-ll 

5.00 


Table 2.3 

Accuracy on ut + u x = 0 with uq(x ) = sin 4 ( 7 rz). 


Method 

N 

Lqo error 

Loo order 

L i error 

L\ order 


20 

1.08e-l 

- 

4.91e-2 

- 


40 

8.90e-3 

3.60 

3.64e-3 

3.75 

WENO-5 

80 

1.80e-3 

2.31 

5.00e-4 

2.86 


160 

1.22e-4 

3.88 

2.17e-5 

4.53 


320 

4.37e-6 

4.80 

6.17e-7 

5.14 


640 

9.79e-8 

5.48 

1.57e-8 

5.30 


20 

5.23e-2 

- 

3.35e-2 

- 


40 

2.47e-3 

4.40 

1.52e-3 

4.46 

CENTRAL-5 

80 

8.32e-5 

4.89 

5.09e-5 

4.90 


160 

2.65e-6 

4.97 

1.60e-6 

4.99 


320 

8.31e-8 

5.00 

4.99e-8 

5.00 


640 

2.60e-9 

5.00 

1.56e-9 

5.00 


In Table 2.2, we show the errors of the fifth order WENO scheme given by the weights (2.58)-(2.59) with the 
smooth indicator (2.63), at time t — 1 for the initial condition u 0 (x) = sin( 7 rx), and compare them with the 
errors of the linear 5 th order upstream central scheme (referred to as CENTRALS in the following tables). 
We can see that fifth order WENO gives the expected order of accuracy starting at about 40 grid points. 

In Table 2.3, we show errors for the initial condition u 0 (a:) = sin 4 ( 7 rx). The order of accuracy for the 
fifth order WENO settles down later than in the previous example. Notice that this is the example for which 
ENO schemes lose their accuracy [63], [67]. 

We emphasize again that the high order conservative finite difference ENO and WENO schemes of third 
or higher order accuracy can only be applied to a uniform grid or a smoothly varying grid, i.e. a grid such 
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that a smooth transformation 


€ = {(*) 

will result in a uniform grid in the new variable £. Here f must contain as many derivatives as the order of 
accuracy of scheme calls for. If this is the case, then (2.65) is transformed to 

ut + £x/(u)f = 0 

and the conservative ENO or WENO derivative approximation is then applied to f{u )%. It is proven in [58] 
that this way the scheme is still conservative, i.e. Lax-Wendroff theorem [51] still applies. 

2.3.3. Boundary conditions. For periodic boundary conditions, or problems with compact support 
for the entire computation (not just the initial data), there is no difficulty in implementing boundary condi- 
tions: one simply set as many ghost points as needed using either the periodicity condition or the compactness 
of the solution. 

Other types of boundary conditions should be handled according to their type: for reflective or symmetry 
boundary conditions, one would set as many ghost points as needed, then use the symmetry/ antisymmetry 
properties to prescribe solution values at those ghost points. For inflow or partially inflow (e.g. a subsonic 
outflow where one of the characteristic waves flows in) boundary conditions, one would usually use the 
physical inflow boundary condition at the exact boundary (for example, if x i is the left boundary and a 
finite volume scheme is used, one would use the given boundary value % as in the monotone flux at 
if xo is the left boundary and a finite difference scheme is used, one would use the given boundary value u *, 
as «o)- Apart from that, the most natural way of treating boundary conditions for the ENO scheme is to use 
only the available values inside the computational domain when choosing the stencil In other words, only 
stencils completely contained inside the computational domain is used in the ENO stencil choosing process 
described in Procedure 2.1. In practical implementation, in order to avoid logical structures to distinguish 
whether a given stencil is completely inside the computational domain, one could set all the ghost values 
outside the computational domain to be very large with large variations (e.g. setting U-j = (10 j) 10 if 
for j = 1,2,..., are ghost points). This way the ENO stencil choosing procedure will automatically 
avoid choosing any stencil containing ghost points. Another way of treating boundary conditions is to use 
extrapolation of suitable order to set the values of the solution in all necessary ghost points. For scalar 
problems this is actually equivalent to the approach of using only the stencils inside the computational 
domain in the ENO procedure. WENO can be handled in a similar fashion. 

Stability analysis (GKS analysis [30], [76]) can be used to study the linear stability when the boundary 
treatment described above is applied to a fixed stencil upwind biased scheme. For most practical situations 
the schemes are linearly stable [3], 

2.3.4. Provable properties in the scalar case. Second order ENO schemes are also TVD (total 
variation diminishing), hence have at least subsequences which converge to weak solutions. There is no 
known convergence result for ENO schemes of degree higher than 2, even for smooth solutions. 

WENO schemes have better convergence results, mainly because their numerical fluxes are smoother. It 
is proven [43] that WENO schemes converge for smooth solutions. Also, Jiang and Yu [44] have obtained an 
existence proof for traveling waves for WENO schemes. This is an important first step towards the proof of 
convergence for shocked cases. 

Even though there are very little theoretical results about ENO or WENO schemes, in practice these 
schemes are very robust and stable. We caution against any attempts to modify the schemes solely for the 
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purpose of stability or convergence proofs. In [69] we gave a remark about a modification of ENO schemes, 
which keeps the formal uniform high order accuracy and makes them stable and convergent for general multi 
dimensional scalar equations. However it was pointed out there that the modification is not computationally 
useful, hence the convergence result has little value. 

The remark in [69] is illustrative hence we reproduce it here. We start with a flux splitting (2.75) 
satisfying (2.76), and notice that the first order monotone scheme 

(2.79) ^ ^ ( f + {ui ) - f + (u i - 1 ) + f-(u i+1 ) - /-(«»)) = 

is convergent (also for multi space dimensions). We now construct a high order ENO approximation in the 
following way: starting from the two point stencil we expand it into a k + 1 point stencil in an 

ENO fashion using the divided differences of / + (it(x)). We then build the A;-th degree polynomial P + (x) 
which interpolates / + (tt(x)) in this stencil. P~{x) is constructed in a similar way, starting from the two 
point stencil {x*, Xi+i}. The scheme is finally defined as 

(2.80) = -£ (P+(x) + P~(x))\ x=xi = R k (u)i 

This scheme is clearly k - th order accurate but is not conservative. We now denote the difference between 
the high order scheme (2.80) and the first order monotone scheme (2.79) by 

(2.81) D(u)i = Rk(u)i - Ri(u)i, 
and limit it by 

(2.82) D{u)i = m(D(u)i , M Ax a ), 

where M > 0 and 0 < a < 1 are constants, and the capping function m is defined by 

a, if \a\ < b ; 
m(a, b) = 6, if a >b\ 

—6, if a < — 6 . 

The modified ENO scheme is then defined by 

du ' - — 

(2.83) ~ = R k {u)i = Ri{u)i + D{u)i. 

We notice that, in smooth regions, the difference between the first order and high order residues, D(u)i , as 
defined in (2.81), is of the size O(Ax), hence the capping (2.82) does not take effect in such regions, if a < 1 
or if a = 1 and M is large enough, when Ax is sufficiently small. This implies that the scheme (2.83) is 
uniformly accurate. Moreover, since 

R k {u) t - Ri{u)i <MAx a 

by (2.82), the high order scheme (2.83) shares every good property of the first order monotone scheme 
(2.79), such as total variation boundedness, entropy conditions, and convergence. From a theoretical point 
of view, this is the strongest result one could possibly hope for a high order scheme. However, the mesh size 
dependent limiting (2.82) renders the scheme highly impractical: the quality of the numerical solution will 
depend strongly on the choice of the parameters M and a, as well as on the mesh size Ax. 
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2.3.5. Systems. We only consider hyperbolic m x m systems, i.e. the Jocobian /'(it) has m real 
eigenvalues 

(2.84) Ai(u) £ ... ^ A m (u) 
and a complete set of independent eigenvectors 

(2.85) n(u), ...,rm(u) . 

We denote the matrix whose columns are eigenvectors (2.85) by 

(2.86) R(u ) = (ri(u), ..., r m (it)) 


Then clearly 

(2.87) R- l {u)f'(u)R(u) = A (it) 

where A(u) is the diagonal matrix with Ai(u), ..., A m (u) on the diagonal. Notice that the rows of i? _1 (it), 
denoted by h(u), ...,Z m (u) (row vectors), are left eigenvectors of f'(u ): 

(2.88) li{u)f'(u) = X i(u)li(u), i = 1, m . 


There are several ways to generalize scalar ENO or WENO schemes to systems. 

The easiest way is to apply the ENO or WENO schemes in a component by component fashion. For the 
finite volume (FV) formulation, this means that we make the reconstruction using ENO or WENO for each 
of the components of u separately. This produces the left and right values at the cell interface 
An exact or approximate Riemann solver, h(u7^ i , u^^), is then used to build the scheme (2.68)-(2.69). 
The exact Riemann solver is given by the exact solution of (2.65) with the following step function as initial 
condition 


(2.89) 


u(x, 0) = 


Vi’ 


i+i 


x < 0; 
x > 0. 


evaluated at the center x = 0. Notice that the solution to (2.65) with the initial condition (2.89) is self- 
similar, that is, it is a function of the variable £ = f , hence is constant along x = 0. If we denote this 
solution by u i+ ^ , then the flux is taken as 


h(u i+h ,« + +i ) = /(« i+ i). 


In the scalar case, the exact Riemann solver gives the Godunov flux (2.70). Exact Riemann solver can be 
obtained for many systems including the Euler equations of compressible gas, which is used very often in 
practice. However, it is usually very costly to get this solution (for Euler equations of compressible gas, an 
iterative procedure is needed to obtain this solution, see [74]). In practice, approximate Riemann solvers 
are usually good enough. As in the scalar case, the quality of the solution is usually very sensitive to the 
choice of approximate Riemann solvers for lower order schemes (first or second order), but this sensitivity 
decreases with an increasing order of accuracy. The simplest approximate Riemann solver (albeit the most 
dissipative) is again the Lax-Friedrichs solver (2.72), except that now the constant a is taken as 


(2.90) 


a = max max I A* (it) I 
u i <j<™ 
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where Aj(u) are the eigenvalues of the Jacobian f f (u ), (2.84). The maximum is again taken over the relevant 
range of u. 

We summarize the procedure in the following 

Procedure 2.6. Component-wise FV ID system ENO and WENO, 

1. For each component of the solution u, apply the scalar ENO Procedure 2.1 or WENO Procedure 
2.2 to reconstruct the corresponding component of the solution at the cell interfaces, , for all i\ 

2. Apply an exact or approximate Riemann solver to compute the flux / I+i for all i in (2.69); 

3. Form the scheme (2.68). 

For the finite difference formulation, a smooth flux splitting (2.75) is again needed. The condition (2.76) 
now becomes that the two Jacobians 

(2 91) d/ + ( M ) df~{u) 

du ’ du 

are still diagonalizable (preferably by the same eigenvectors R(u) as for and have only non-negative 

/ non-positive eigenvalues, respectively. We again recommend the Lax-Friedrichs flux splitting (2.77), with 
a given by (2.90), because of its simplicity and smoothness. A somewhat more complicated Lax- Friedrichs 
type flux splitting is: 

/*( u ) = ± #(«) A R~ x (u)u), 

where R(u) and J?“ 1 (u) are defined in (2.86), and 

A = diagi Ai,...,A m ) 

where A^ = max u |Aj(k)|, and the maximum is again taken over the relevant range of u. This way the 
dissipation is added in each field according to the maximum size of eigenvalues in that field, not globally. 
One could also use other flux splittings, such as the van Leer splitting for gas dynamics [79]. However, for 
higher order schemes, the flux splitting must be sufficiently smooth in order to retain the order of accuracy. 

With these flux splittings, we can again use the scalar recipes to form the finite difference scheme: just 
compute the positive and negative fluxes i and f~ i component by component. 

We summarize the procedure in the following 

Procedure 2.7. Component-wise FD ID system ENO and WENO. 

1. Find a flux splitting (2.75). The simplest example is the Lax-Friedrichs flux splitting (2.77), with a 
given by (2.90); 

2. For each component of the solution u , apply the scalar Procedure 2.5 to reconstruct the corresponding 
component of the numerical flux / i+ ± ; 

3. Form the scheme (2.73). 

These component by component versions of ENO and WENO schemes are simple and cost effective. 
They work reasonably well for many problems, especially when the order of accuracy is not high (second or 
sometimes third order). However, for more demanding test problems, or when the order of accuracy is high, 
we would need the following more costly, but much more robust characteristic decompositions. 

To explain the characteristic decomposition, we start with a simple example where f(u) = An in (2.65) 
is linear and A is a constant matrix. In this situation, the eigenvalues (2.84), the eigenvectors (2.85), and 
the related matrices R, R 1 and A (2.86)- (2.87), are all constant matrices. If we define a change of variable 

(2.92) v — R~ l 
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then the PDE (2.65) becomes diagonal: 

(2.93) v t + Av* = 0 

that is, the m equations in (2.93) are decoupled and each one is a scalar linear convection equation of the 
form 


(2.94) w t + AjtUx = 0. 

We can thus use the reconstruction or flux evaluation techniques for the scalar equations, discussed in 
Sections 2.3.1 and 2.3.2, to handle each of the equations in (2.94). After we obtain the results, we can “come 
back” to the physical space u by using the inverse of (2.92): 


(2.95) 


u = Rv 


For example, if the reconstructed polynomial for each component j in (2.93) is denoted by Qj{x), then we 
form 


( 9i (*) \ 

(2.96) <l( x ) = 

\ Qm{x) ) 

and obtain the reconstruction in the physical space by using (2.95): 

(2.97) p ( x ) = Rq{x) 


The flux evaluations for the finite difference schemes can be handled similarly. 

We now come to the situation where f(u) is not constant. The trouble is that now all the matrices R{u), 
R~ l (u) A(u) are dependent upon u. We must “freeze” them locally in order to carry out a similar procedure 
as in the constant coefficient case. Thus, to compute the flux at the cell boundary we would need an 

approximation to the Jocobian at the middle value u i+ i. This can be simply taken as the arithmetic mean 

( 2 . 98 ) u i+ i = i (Ui + tti+i) , 

or as a more elaborate average satisfying some nice properties, e.g. the mean value theorem 


(2.99) f(u i+ i)~ f(ui) = f'(u i+ i)(u i+ i-ui). 

Roe average [62] is such an example for the compressible Euler equations of gas dynamics and some other 
physical systems. It is also possible to use two different one-sided Jacobians at a higher computational cost 

[25]. 

Once we have this u i+i , we will use R(u t+ i), /T l (u i+i ) and A(u i+i ) to help evaluating the numerical 
flux at i <+ i. We thus omit the notation i + ^ and still denote these matrices by R, R 1 and A, etc. We 
thpn repeat the procedure described above for linear systems. The difference here being, the matrices R, 
r-i and A are different at different locations x t+ i, hence the cost of the operation is greatly increased. 
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In summary, we have the following procedures: 

Procedure 2.8. Characteristic- wise FV ID ENO and WENO. 

1. Compute the divided or undivided differences of the cell averages tZ, for all i ; 

2. At each fixed a^+i, do the following: 

(a) Compute an average state , using either the simple mean (2.98) or a Roe average satisfying 
(2.99); 

(b) Compute the right eigenvectors, the left eigenvectors, and the eigenvalues of the Jacobian 
/'(Uj+i ), (2.84)- (2.87), and denote them by 

R = R(u i+i ), i), A = A( Ui+i ); 

(c) Transform all those differences computed in Step 1, which are in the potential stencil of the 
ENO and WENO reconstructions for obtaining uj^i, to the local characteristic fields by using 
(2.92). For example, 

Vj = R~ l uj, j in a neighborhood of i; 

(d) Perform the scalar ENO or WENO reconstruction Procedure 2.3, for each component of the 
characteristic variables u, to obtain the corresponding component of the reconstruction 

(e) Transform back into physical space by using (2.95): 




3. Apply an exact or approximate Riemann solver to compute the flux / i+i for all i in (2.69); then 
form the scheme (2.68). 

Similarly, the procedure to obtain a finite difference ENO-Roe type scheme using the local characteristic 
variables is: 


Procedure 2.9. Characteristic-w ise FD ID system, Roe-type. 

1. Compute the divided or undivided differences of the flux f(u) for all i\ 

2. At each fixed do the following: 

(a) Compute an average state , using either the simple mean (2.98) or a Roe average satisfying 
(2.99); 

(b) Compute the right eigenvectors, the left eigenvectors, and the eigenvalues of the Jacobian 
f r {u i+ i) y (2.84)-(2.87), and denote them by 

R — R(u i+ i) y R =R (u i+ i), A = A(u i+ i); 

(c) Transform all those differences computed in Step 1, which are in the potential stencil of the 
ENO and WENO reconstructions for obtaining the flux / i+ i , to the local characteristic fields 
by using (2.92). For example, 

Vj — R~ l f(uj ), j in a neighborhood of i; 

(d) Perform the scalar ENO or WENO Roe-type Procedure 2.4, for each component of the char- 
acteristic variables u, to obtain the corresponding component of the flux v i+ i. The Roe speed 

1 is replaced by the eigenvalue Aj(u i+ i ) for the l-th. component of the characteristic variables 
v\ 
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(e) Transform back into physical space by using (2.95): 

£+§=**«+* 

3. Form the scheme (2.73). 

Finally, the procedure to obtain a finite difference flux splitting ENO or WENO scheme using the local 
characteristic variables is: 


Procedure 2.10. Characteristic- wise FD ID system, flux splitting. 

1. Compute the divided or undivided differences of the flux f(u) and the solution u for all i; 

2. At each fixed i, do the following: 

' 2 _ 

(a) Compute an average state u i+ ^ , using either the simple mean (2.98) or a Roe average satisfying 
(2.99); 

(b) Compute the right eigenvectors, the left eigenvectors, and the eigenvalues of the Jacobian 
f'(u i+ i), (2.84)-(2.87), and denote them by 

R = R(u i+ ±), R~ 1 = R~ 1 {u i+ i), A = A(u i+ i); 

(c) Transform all those differences computed in Step 1, which are in the potential stencil of the 
ENO and WENO reconstructions for obtaining the flux f i+ i , to the local characteristic fields 
by using (2.92). For example, 

i; j =R _1 u J , g j = R~ l f{u j ) , j in a neighborhood of i; 

(d) Perform the scalar flux splitting ENO or WENO Procedure 2.5, for each component of the 
characteristic variables, to obtain the corresponding component of the flux gf + ±. For the 
most commonly used Lax-Friedrichs flux splitting, we can use, for the 1-th component of the 
characteristic variables, the viscosity coefficient 


a = max IM“j)ll 

1<J<N 

(e) Transform back into physical space by using (2.95): 




3. Form the flux by taking 

•^*+4 ^*+4 ^ 

and then form the scheme (2.73). 

There are attempts recently to simplify this characteristic decomposition. For example, for the com- 
pressible Euler equations of gas dynamics, Jiang and Shu [43] used smooth indicators based on density and 
pressure to perform the so-called pseudo characteristic decompositions. There are also second and sometimes 
third order component ENO type schemes [56], [54], with limited success for higher order methods. 
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3. Multi Space Dimensions. 

3.1. Reconstruction and Approximation in Multi Dimensions. In this section we describe how 
the ideas of reconstruction and approximation in Sect. 2.1 are generalized to multi space dimensions. We 
will concentrate our discussion in 2D, although things carry over to higher dimensions as well. 

In most part of this section we will consider Cartesian grids, that is, the domain is a rectangle 

(3-1) [a, b] x [c, d\ 

covered by cells 

( 3 - 2 ) = [ari_i,x i+ i] x [Vj_i,y j+ ^ 1 < i < N x , 1 < j < N y 

where 

a = X \ < X % < < x iV a: + ^ = ^ 


and 


c = v\ <»§<■■•< vn v ~\ < y Nv+ 1 = d. 

The centers of the cells are 

(3.3) 

and we still use 

(3.4) 
and 

( 3 ‘ 5 ) A Vj =y j+h - yj _ h , j = 1, 2, ..., N y 

to denote the grid sizes. We denote the maximum grid sizes by 


fo, Vj), x t = - (*<_J + * i+ i) , Vj = ~ {y 3 _^ + y j+ ! ) , 


Ax z — % ~ 1» 2? N x 


(3.6) 


Ax = max Ax 4, Ay = max Ay,, 

1 <i<N* * 1 <j<N v 


and assume that Ax and Ay are of the same magnitude (their ratio is bounded from above and below during 
refinement). Finally, 


(3.7) 


A = max (Ax, Ay). 


3.1.1. Reconstruction from cell averages. The approximation problem we will face, in solving 
hyperbolic conservation laws using cell averages (finite volume schemes, see Sect. 3.3), is still the following 
reconstruction problem. 

Problem 3.1. Two dimensional reconstruction. 

Given the cell averages of a function v(x,y): 


(3.8) 


i r y i + i 


i = 1,2,...,N X , j = 1, 2, ..., N y , 
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find a polynomial Pij(x,y), preferably of degree at most k - 1, for each cell I <j, such that it is a k - th order 
accurate approximation to the function t>(x, y) inside I if. 

(3.9) Pijfay) = v{x,y)+0( A fc ), (x,y) € i = 1, j = 1, JV*. 

In particular, this gives approximations to the function u(x,y) at the cell boundaries 



= Pa( x i+^y)’ v t-\, y = i = 1 > •••’ 

N x , <y< y j+ 1 


v~ j+i = Pij(x, = Pij(x, = 1. 

<x<x i+ l 

which are 

fc-th order accurate: 


(3.10) 

v t+i,y = V ( x i+^y) +0(A k ), i -- 0,1 ,...,N X , 

+ 

d? 

VI 

VI 

and 



(3.11) 

V Xj+\ =V ( X >yj+$) + °( Ak ')' j = 0> 1) ■"! Ny, 

X i-\ — X — X i+\ * 


Again we will not discuss boundary conditions in this section. We thus assume that Vij is also available 
for i < 0, i > N x and for j < 0, j > N y if needed. 

In the following we describe the general procedure to solve Problem 3.1. 

Given the location I tj and the order of accuracy k y we again first choose a “stencil”, based on 
neighboring cells, the collection of these cells still being denoted by We then try to find a polynomial 

of degree at most k - 1, denoted by p(x, y) (we again drop the subscript ij when it does not cause confusion), 
whose cell average in each of the cells in S(z, j) agrees with that of v(x, y ): 

l f V rn+% f X l+^ _ 

(3.12) -T — -T / / p{Z,v)d£drj = v lm , if hm eS(i,j ). 

Ax f Ay m Jy^ ^ J Xi _^ 

We first remark that there are now many more candidate stencils S(i, j) than in the ID case, More 
importantly, unlike in the ID case, here we encounter the following essential difficulties: 

• Not all of the candidate stencils can be used to obtain a polynomial p(x, y) of degree at most k- 1 
satisfying condition (3.12). 

For example, it is an easy exercise to show that neither existence nor uniqueness holds, if one wants 
to reconstruct a first degree polynomial p(x, y) satisfying (3.12) for the three horizontal cells 

= {Ii-ij, I%j , /i-j-ij} • 

To see this, let’s assume that 

h-U = [-2A, - A] x [0, A] ,Ij j = [-A,0] x [0, A], 
li+i,} = [0, A] x [0, A], 

and the first degree polynomial p(x, y) is given by 

p(x, y) = a + fix + 73/ 

then condition (3.12) implies 


a - %A(3 + | A7 — 


a — gA j3 + 5A7 = 

Vi } j 

a+\&.(3+\ki = 



which is a singular linear system for a, ft and 7. 
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• Even if one obtains such a polynomial p(x , y), there is no guarantee that the accuracy conditions 
(3.9) will hold. We again use the same simple example. If we pick the function 

v{x,y) = 0 , 

then one of the polynomials of degree one satisfying the condition (3.12) is 

p(x, y) = A - 2y 

clearly the difference 


v(x, 0) - p(x , 0) = — A 

is not at the size of 0( A 2 ) in < x < , as is required by (3.9). 

This difficulty will be more profound for unstructured meshes such as triangles. See, for example, [l]. 
For rectangular meshes, if we use the tensor products of ID polynomials, i.e. use polynomials in Q k ~ 1 : 

k - 1 fc-l 

p(x, y) = EE QlmX y 

m =0 1=0 

then things can proceed as in ID. We restrict ourselves in the following tensor product stencils: 


S r8 (h j) — {hm -i — r<l<i + k — 1— r, j — s<m<j + k- l-s} 


then we can address Problem 3.1 by introducing the two dimensional primitives: 


Clearly 


V(x,y) 



v(€,v)<%dri. 


f y j+i 

V(x i+ i,y i+ i)= / v(£, r])d£dr) 

J — oc J — oo 

j j 

= v lrn Ax[Ay m , 

771= — OO 1 = — OO 

hence as in the ID case, with the knowledge of the cell averages v we know the primitive function V exactly 
at cell corners. 

On a tensor product stencil 


Srs(i,j) = {(*<+$, Vm+l) :*-r-l<Z<* + fc-l-r, j-s-l<m<j + k~l-s} 


there is a unique polynomial P(x,y) in Q k which interpolates V at every point in S r8 (i,j). We take the 
mixed derivative of the polynomial P to get: 


p(x,y) 


d 2 P(x,y) 

dxdy 


then p(x , y) is in Q k approximates v(x , y ), which is the mixed derivative of V(x, y), to k-th order: 


v(x,y) -p(x,y) = 0{A k ) 
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and also satisfies (3.12): 


1 _ — f m+ i f ‘ +i dr) 

Ax l AymJ Vm _ i J x ,_ h 


1 f Vm+ i f 1 ' 

iAy m Jy m _i ^ T i- 


i-i &P_ 

d^dr] 


(£,»? )d£d.T) 


i •"<-* 

(P(x l+ ±,y m+ 1 ) - F(x l+ i,y m _i) 


Ax/ Ay 

_ 1 
AxiAy m 

-P(x,_i,y m+ i) + P(x,_i,y m _i)) 

= ■ nxi+ * ,y ’"-* ) 

- V(®, _ i , y m+ 1 ) + V (x,_ i , y m _ i )) 

= -T— )r — f m+i f 1+i v(t,T) )d£dr) = vim , 

Ax/Ay m Jy m _ ^ Jx l _^ 

i-r <1 <i + k-\-r, j-s<m<j + k-~l-s. 


This gives us a practical way to perform the reconstruction in 2D. We first perform a one dimensional 
reconstruction (Problem 2.1), say in the y direction, obtaining one dimensional cell averages of the function 
v in the other direction (say in the x direction). We then perform a reconstruction in the other direction. 

It should be remarked that the cost to do this 2D reconstruction is very high: For each grid point, if 
the cost to perform a one dimensional reconstruction is c, then we need 2c per grid point to perform this 2D 
reconstruction. In general n space dimensions, the cost grows to nc. 

We also remark that to use polynomials in Q h ~ l is a waste: to get the correct order of accuracy only 
polynomials in P k ~ l is needed. However, there is no natural way of utilizing polynomials in P k 1 (see the 
comments above and also the paper of Abgrall [1]). 

The reconstruction problem, Problem 3.1, can also be raised for general, non-Cartesian meshes, such as 
triangles. However, the solution becomes much more complicated. For discussions, see for example [1]. 

3.1.2. Conservative approximation to the derivative from point values. The second approx- 
imation problem we will face, in solving hyperbolic conservation laws using point values (finite difference 
schemes, see Sect. 3.2), is again the following problem in obtaining high order conservative approximation 
to the derivative from point values [69, 70]. As in the ID case, here we also assume that the grid is uniform 
in each direction. We again ignore the boundary conditions and assume that Vij is available for i < 0 and 
i > N x , and for j < 0 and j > N y . 


Problem 3.2. Two dimensional conservative approximation. 

Given the point values of a function v(x, y): 

(3.13) Vij = v{xi,yj), i = 1,2,..., N xy j = 1,2, 
find numerical flux functions 

(3.14) 1 j = > *•■? * — 1,...,AT X 

and 

(3.15) Vi,j+± =v(v itj - t> ..., v it j+ k-i-s), j = 0, 1, •••, N y 
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such that the flux differences approximate the derivatives v x (x, y) and v y (x, y) to k-th order accuracy: 


(3.16) 

Ax i Vi+ i’ j v, -bi) 

= v x (xi,yj) + 0(Ax k ), 

4 

i— i 

cT 

II 

and 




(3.17) 

Ay ~ Vi ’3~i) 

= v y{ x i, Vj) + 0(Ay k ), 

a? 

r—4 

o 

II 


The solution of this problem is essential for the high order conservative schemes based on point values 
(finite difference) rather than on cell averages (finite volume). 

Having seen the complication of reconstructions in the previous subsection for multi space dimensions, 
it is a good relieve to see that conservative approximation to the derivative from point values is as simple in 
multi dimensions as in ID. In fact, for fixed j , if we take 

w{x) = v(x,yj) 

then to obtain v x (%i,yj) = w f (xi) we only need to perform the one dimensional procedure in Sect. 2.1, 
Problem 2.2, to the one dimensional function w(x). Same thing for v y (x, y). 

As in the ID case, the conservative approximation to derivatives, of third order accuracy or higher, can 
only be applied to uniform or smoothly varying meshes (curvilinear coordinates). It cannot be applied to 
general unstructured meshes such as triangles, unless conservative is given up. 

3.2. ENO and WENO Approximations in Multi Dimensions. For solving hyperbolic conserva- 
tion laws in multi space dimensions, we are again interested in the class of piecewise smooth functions. We 
define a piecewise smooth function v(x, y) to be such that, for each fixed y, the one dimensional function 
w(x) = v(x,y) is piecewise smooth in the sense described in Sect. 2.2. Likewise, for each fixed x , the 
one dimensional function w(y) = v(x , y ) is also assumed to be piecewise smooth. Such functions are again 
“generic” for solutions to multi dimensional hyperbolic conservation laws in practice. 

In the previous section, we have already discussed the problems of reconstruction and conservative 
approximations to derivatives in multi space dimensions. At least for the Cartesian type grids, both the 
reconstruction and the conservative approximation can be obtained from one dimensional procedures. 

For the reconstruction, we first use a one dimensional ENO or WENO reconstruction, Procedure 2.1 
or 2.2, on the two dimensional cell averages, say in the y direction, to obtain one dimensional cell averages 
in x only. Then, another one dimensional reconstruction in the remaining direction, say in the x direction, 
is performed to recover the function itself, again using the one dimensional ENO or WENO methodology, 
Procedure 2.1 or 2.2. 

For the conservative approximation to derivatives, since they are already formulated in a dimension by 
dimension fashion, one dimensional ENO and WENO procedures can be trivially applied. In effect, the 
FORTRAN program for the 2D problem is the same as the one for the ID problem, with an outside “do 
loop”. 

What happens to general geometry which cannot be covered by a Cartesian grid? 

If the domain is smooth enough, it usually can be mapped smoothly to a rectangle (or at least to a union 
of non-overlapping rectangles). That is, the transformation 

( 318 ) Z = €(x,y), T] = r)(x,y) 
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maps the physical domain D where (x, y) belongs, to a rectangular computational domain 

(3.19) a < £ < 6, c<r]<d. 

We require the transformation functions (3.18) to be smooth (i.e. it has as many derivatives as the accuracy 
of the scheme calls for). Using chain rule, we could write, for example, 

(3.20) v x = £ x vs + VxV v 

We can then use our ENO or WENO approximations on and v v , as they are now defined in rectangular 
domains. The smoothness of £ x and 7] x will guarantee that this leads to a high order approximation to v x 
as well through (3.20). 

If the domain is really ugly, or if one wants to use unstructured meshes for other purposes (e.g. for 
adaptivity), then ENO and WENO approximations for unstructured meshes must be studied. This is a 
much less matured subject at present. We refer the readers to [1] for some efforts in this direction. 

3.3* ENO and WENO Schemes for Multi Dimensional Conservation Laws. In this section we 
describe the ENO and WENO schemes for 2D conservation laws: 


(3.21) 


u t (x,y,t) + fx(u(x,y,t)) + g y {u(x,y,t)) = 0 


again equipped with suitable initial and boundary conditions. 

Although we present everything in 2D, most of the discussion is also valid for higher dimensions. 

We again concentrate on the discussion of spatial discretizations, and will leave the time variable t 
continuous (the method- of-lines approach). Time discretizations will be discussed in Sect. 4.2. 

In most of the discussion in this section, our computational domain is rectangular, given by (3.1). Our 
grids will thus be Cartesian, given by (3.2) and (3.3). Unstructured meshes will only be mentioned briefly. 

We do not discuss boundary conditions in this section. We thus assume that the values of the numerical 
solution are also available outside the computational domain whenever they are needed. This would be the 
case for periodic or compactly supported problems. Two dimensional boundary condition treatments are 
similar to the one dimensional case discussed in Sect. 2.3.3. 


3.3.1. Finite volume formulation in the scalar case. For finite volume schemes, or schemes based 
on cell . averages, we do not solve (3.21) directly, but its integrated version. We integrate (3.21) over the 
interval in to obtain 


(3.22) 


where 

(3.23) 



is the cell average. We approximate (3.22) by the following conservative scheme 


(3.24) 


duij (t) 

dt 


1 ft f 

\ 1 /. . \ 

~Kxi v i+ ^ Ji ~^. 

) A Vj 9l 'i~l) 
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where the numerical flux f i+ ij is defined by 


^ 3 ' 25 ) fi+l’i ^2 w ah { U i+± t y j+ p a Ay i ' U t+$ ,» J +/J e ,A» i ) > 

a 

where /3 a and w a are Gaussian quadrature nodes and weights, for approximating the integration in y: 

1 f y j+i 

— J f(u(x i+ L,y,t)dy 


inside the integral form of the PDE (3.22), and 1 are the A;~th order accurate reconstructed values 

1 i 2 yV 

obtained by ENO or WENO reconstruction described in the previous section. As before, the superscripts db 
imply the values are obtained within the cell I{j (for the superscript -) and the cell Ii+ij (for the superscript 
-fi), respectively. The flux 9i t3+ i is defined similarly by 


( 3 - 26 ) ki+h 

a 

for approximating the integration in x: 


+p a AxiJ+±’ U Xi+/3 Q Ax, 


ij+i ) ’ 


1 f X i + 1 

— / 2 g(u{x,y j+ ±,t)dx 

* x ' •/«,_$ “ * 

inside the integral form of the PDE (3.22). tt* A are again the k-th order accurate reconstructed values 
obtained by ENO or WENO reconstruction described in the previous section, h is again the one dimensional 
monotone flux, examples being given in (2. 70)- (2. 72). 

We summarize the procedure to build a finite volume ENO or WENO 2D scheme (3.24), given the cell 
averages {u i3 } (we again drop the explicit reference to the time variable t), and a one dimensional monotone 
flux h ) as follows: 


Procedure 3.1. Finite volume 2D scalar ENO and WENO, 

1. Follow the procedures described in Sect. 3.2, to obtain ENO or WENO reconstructed values at 
the Gaussian points, u f + i iyj +p aAy an d U x i +p a Ax i j+^’ ^°^ ce that this step involves two one 
dimensional reconstructions, each one to remove a one dimensional cell average in one of the two 
directions. Also notice that the optimal weights used in the WENO reconstruction procedure are 
different for different Gaussian points indexed by a; 

2. Compute the flux /i+i j and 9ij+i using (3.25) and (3.26); 

3. Form the scheme (3.24). 


We remark that the finite volume scheme in 2D, as described above, is very expensive due to the following 
reasons: 

• A two dimensional reconstruction, at the cost of two one dimensional reconstruction per grid point, 
is needed. For general n space dimensions, the cost becomes ti one dimensional reconstruction per 
grid point; 

• More than one quadrature points are needed in formulating the flux (3.25)-(3.26), for order of 
accuracy higher than two. Thus, for ENO, although the stencil choosing process needs to be done 
only once, the reconstruction (2.10) has to be done for each quadrature point used in the flux 
formulation. For WENO, the optimal weights are also different for each quadrature point. This 
becomes much more costly for n > 2 dimension, as then the fluxes are defined by integrals in n - 1 
dimension and a n — 1 dimensional quadrature rule must be used. 
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This is why multidimensional finite volume schemes of order of accuracy higher than 2 are rarely used. 
For 2D, based on [34], Casper [11] has coded up a fourth order finite volume ENO scheme for Cartesian 
grids, sec also [12]. 3D finite volume ENO code of order of accuracy higher than 2 does not exist yet, to the 
author’s knowledge. 

At the second order level, the cost is greatly reduced because: 

• There is no need to perform a reconstruction, as the cell average U{j agrees with the point value at 

the center u(xi,yj) to second order 0(A 2 ); 

• The quadrature rule in defining the flux (3.25)-(3.26) needs only one (mid) point. 

One advantage of finite volume ENO or WENO schemes is that they can be defined on arbitrary meshes, 
provided that an ENO or WENO reconstruction on that mesh is available. See, for example, [1]. 

3.3.2. Finite difference formulation in the scalar case. Here we assume a uniform grid and solve 
(3.21) directly using a conservative approximation to the spatial derivative: 


(3.27) 


diiij (£) 

dt 


1 

Ax 




1 

Ay 



“ 9ij- 


1 

2 


) 


where u x] (<) is the numerical approximation to the point value u(xi,yj,t). 

The numerical flux / l+i j is obtained by the one dimensional ENO or WENO approximation procedure, 
Procedure 2.4 or 2.5, with v(x) = f(u(x,yj,t)) and with j fixed. Likewise, the numerical flux g ij+ i is 
obtained by the one dimensional ENO or WENO approximation procedure, with v(y) = f{u(x t ,y,t)) and 
with i fixed. 

All the one dimensional discussions in Sect. 2.3.2, such as upwinding, ENO-Roe, flux splitting, etc., can 
be applied here dimension by dimension. 

The discussion here is also valid for higher spatial dimension n. In effect, it is the same one dimensional 
conservative derivative approximation applied to each space dimension. 

It is a straight forward exercise [13] to show that, in terms of operation count, the finite difference ENO 
or WENO schemes are about a factor of 4 less than the finite volume counterpart of the same order. In 3D 
this factor becomes about 9. 

We thus strongly recommend the usage of the finite difference version of ENO and WENO schemes (also 
called ENO and WENO schemes based on point values), whenever possible. 

3.3.3. Provable properties in the scalar case. Second order ENO schemes are also maximum norm 
non-increasing. Of course, this stability is too weak to imply any convergence. As was mentioned before, 
there is no known convergence result for ENO schemes of order higher than 2, even for smooth solutions. 

WENO schemes have better convergence results also in the current multi-D case, mainly because their 
numerical fluxes are smoother. It is proven [43] that WENO schemes converge for smooth solutions. 

We again emphasize that, even though there are very little theoretical results about ENO or WENO 
schemes, in practice they are very robust and stable. We once again caution against any attempts to modify 
the schemes solely for the purpose of stability or convergence proofs. In fact the modification of ENO 
schemes in [69], presented in Sect. 2.3.4, which keeps the formal uniform high order accuracy, actually 
produces schemes which arc convergent to entropy solutions for general multi dimensional scalar equations. 
However it was pointed out there that the modification is not computationally useful, hence the convergence 
result has little value. 

3.3.4. Systems. The advice here is that, when the fluxes are computed along a cell boundary, a one 
dimensional local characteristic decomposition normal to the boundary is performed. Also, the monotone 
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flux is replaced with a one dimensional exact or approximate Riemann solver. Thus, the discussion in Sect. 
2.3.5 can be applied here. 

There are discussions in the literature about truly multi-dimensional recipes. However, these tend to 
become extremely complicated for order of accuracy higher than two, so they have not been used extensively 
in practice. Another reason to suggest against using such complicated truly multidimensional recipes for 
order of accuracy higher than two is that, while dimension by dimension schemes as advocated in these 
lecture notes are not rotationally invariant, the direction related non-symmetry actually diminishes with 
increased order [13]. 

4. Further Topics. 


4.1. Further Topics in ENO and WENO Schemes. In this section we discuss some miscellaneous 
(but not necessarily unimportant!) topics in ENO and WENO schemes. 


4.1.1. Subcell resolution. This idea was first raised by Harten [35]. The observation is that, since 
in interpolating the primitive V, two points must be included in the initial stencil (see Procedure 2.1), one 
cannot avoid having at least one cell for each discontinuity, inside which the reconstructed polynomial is 
not accurate (0(1) error there). We can clearly see this 0(1) error in the ENO interpolation in Figure 2.1. 
The reconstruction in this shocked cell, although inaccurate, will always be monotone (Property 2 in Sect. 
2.2.1), so stability will not be a problem. However, it does cause a smearing of the discontinuity (over one 
cell, initially). 

If we are solving a truly nonlinear shock, then characteristics flow into the shock, thus any error one 
makes during time evolution tends to be absorbed into the shock (we also say that the shock has a self 
sharpening mechanism). However, we are less lucky with a linear discontinuity, such as a discontinuity 
carried by the linear equation ut + u x = 0. Such linear discontinuities are also called contact discontinuities 
in gas dynamics. The characteristics for such cases are parallel to the discontinuity, hence any numerical 
smearing tends to accumulate and the discontinuity becomes progressively more smeared with time (Harten 
argues that the smearing of the discontinuity is at the rate of 0{ Ax 1- ^ 1 ) where k is the order of the scheme. 
Although higher order schemes have less smearing, when time is large the smearing is still very significant. 

Harten [35] makes the following simple observation: in the shocked cell instead of using the reconstruc- 
tion polynomial Pi{x), which is highly inaccurate (the only useful information it carries is the cell average 
in the cell), one could try to find the location of the discontinuity inside the cell /*, say at x s , and then use 
the neighboring reconstructions Pi-\{x) extended to x s from left and Pt+i(x) extended to x 3 from right. To 
find the shock location, one could argue that pi_ i(x) is a very accurate approximation to v{x) up to the 
discontinuity x s from left, and Pi+\{x) is a very accurate approximation to u(ar) up to the discontinuity x 3 
from right. We thus extend Pi-\{x) from the left into the cell /*, and extend Pi + i(x) from the right into the 
cell Jj, and require that the cell average be preserved: 


(4.1) 



Pi+i{x)dx 


A XiVi. 


It can be proven that under very general conditions, (4.1) has only one root x 8 inside the cell /», hence one 
could use Newton iterations to find this root. 

Subcell resolution can be applied to both finite volume and finite difference ENO and WENO schemes 
[35], [70], However, it should be applied only to sharpen contact discontinuities. It is quite dangerous to 
apply subcell resolution to a shock, since it might generate entropy violating expansion shocks in the 
numerical solution. 
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Another very serious restriction about subcell resolution is that it is very difficult to be applied to 2D. 
However, see Siddiqi, Kimia and Shu [73], where a geometrical ENO is used to extend the subcell resolution 
idea to 2D for image processing problems (we termed it geometric ENO, or GENO). 

4.1.2. Artificial compression. Another very useful idea to sharpen a contact discontinuity is the 
artificial compression, first developed by Harten [32] and further improved by Yang [83]. The idea is to 
increase the magnitude of the slope of a reconstruction, of course subject to certain monotonicity restrictions, 
near such a discontinuity. Notice that this goes against the idea of limiting, which typically decreases the 
magnitude of the slope of a reconstruction. 

Artificial compression can be applied both to finite volume and to finite difference ENO and WENO 
schemes [83], [70], [43]. Unlike subcell resolution, artificial compression can also be applied easily to multi 
space dimensions, at least in principle. 

4.1.3. Other building blocks. It is not necessary to stay within polynomial building blocks, although 
polynomials arc the most natural functions to work with. For some applications, other building blocks, such 
as rational functions, trigonometric polynomials, exponential functions, radial functions, etc., may be more 
appropriate. The idea of ENO or WENO can be applied also in such situations. The key idea is to find 
suitable “smooth indicators”, similar to the Newton divided differences for the polynomial case, for applying 
the ENO or WENO idea. See [16] and [40] for some examples. 

4.2. Time Discretization. Up to now we have only considered spatial discretizations, leaving the 
time variable continuous (method of lines). In this section we consider the issue of time discretization. The 
techniques discussed in this section can also be applied to other types of spatial discretizations using the 
method of lines approach, such as various TVD and TVB schemes [52, 78, 65] and discontinuous Galerkin 
methods [18, 19, 20, 21, 17]. 

4.2.1. TVD Runge-Kutta methods. A class of TVD (total variation diminishing) high order Runge- 

Kutta methods is developed in [69] and further in [29] . 

These Runge-Kutta methods are used to solve a system of initial value problems of ODEs written as: 

(4.2) u t = L ( u )> 

resulting from a method of lines spatial approximation to a PDE such as: 

(4.3) 

We have written the equation in (4.3) as a ID conservation law, but the discussion which follows apply to 
general initial value problems of PDEs in any spatial dimensions. Clearly, L{ u) in (4.2) is an approximation 
(e.g. ENO or WENO approximation in these lecture notes), to the derivative -f(u) x in the PDE (4.3). 

If we assume that a first order Euler forward time stepping: 

(4.4) u n + l =u n + AtL(u n ) 
is stable in a certain norm: 

(4.5) ll« n+1 ll < ll« n ll 
under a suitable restriction on At: 

(4.6) At<Ati, 


42 



then we look for higher order in time Runge-Kutta methods such that the same stability result (4.5) holds 
under a perhaps different restriction on At: 


(4-7) At<cAti. 

where c is termed the CFL coefficient for the high order time discretization. 

We remark that the stability condition (4.5) for the first order Euler forward in time (4.4) is easy to 
obtain in many cases, such as various TVD and TVB schemes in ID (where the norm is the total variation 
norm) and in multi dimensions (where the norm is the L°° norm), see, e.g. [52, 78, 65], 

Originally in [69, 66] the norm in (4.5) was chosen to be the total variation norm, hence the terminology 
“TVD time discretization”. 

As it stands, the TVD high order time discretization defined above maintains stability in whatever norm, 
of the Euler forward first order time stepping, for the high order time discretization, under the time step 
restriction (4.7). For example, if it is used for multi dimensional scalar conservation laws, for which TVD 
is not possible but maximum norm stability can be maintained for high order spatial discretizations plus 
forward Euler time stepping (e.g. [20]), then the same maximum norm stability can be maintained if TVD 
high order time discretization is used. As another example, if an entropy inequality can be proved for the 
Euler forward, then the same entropy inequality is valid under a high order TVD time discretization. 

In [69], a general Runge-Kutta method for (4.2) is written in the form: 

( 4 - 8 ) « (i) = J2 (mku w + AtAkL(u<*>)) , i= 1 , ..., m 

k = o v 


Clearly, if all the coefficients are nonnegative a ik > 0, 0 ik > 0, then (4.8) is just a convex combination of 
the Euler forward operators, with At replaced by J^Af, since by consistency a lk = 1. We thus have 

Lemma 4,1. [69] The Runge-Kutta method (4.8) is TVD under the CFL coefficient (4.7): 

(4-9) c = min^p^, 

».* 0ih 

provided that a**, > 0, (3^ > 0. 


In [69], schemes up to third order were found to satisfy the conditions in Lemma 4.1 with CFL coefficient 
equal to 1. 

The optimal second order TVD Runge-Kutta method is given by [69, 29]: 

(4.10) w (!) = u n + A tL(u n ) 

<*" +1 = ]«“ + j« (1) + 1 

with a CFL coefficient c = 1 in (4.9). 

The optimal third order TVD Runge-Kutta method is given by [69, 29]: 

u (i) — u n + A tL(u n ) 

(4-11) u (2 > = + 1 U (1) + 1 

u n+1 = \u n + + ?A tL( u W), 
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with a CFL coefficient c — 1 in (4.9). 

Unfortunately, it is proven in [29] that no four stage, fourth order TVD Runge-Kutta method exists with 
nonnegative a ik and (3 lk . We thus have to consider the situation where a ik > 0 but Afc might be negative. 
In such situations we need to introduce an adjoint operator L. The requirement for L is that it approximates 
the same spatial derivative(s) as L, but is TVD (or stable in another relevant norm) for first order Euler, 
backward in time: 

(4.12) « n+1 = «" " 

This can be achieved, for hyperbolic conservation laws, by solving the backward in time version of (4.3): 

(4.13) u t = /( u )x- 

Numerically, the only difference is the change of upwind direction. Clearly, L can be computed with the 
same cost as that of computing L. We then have the following lemma: 


Lemma 4.2. [69] The Runge-Kutta method (4.8) is TVD under the CFL coefficient (4.7): 


(4.14) 


c 


= min 
i,k 


otik 

lAfcl’ 


provided that cti k > 0, and L is replaced by L for negative /?*fc - 


Notice that, if for the same fc, both L(u «) and L(u^) must be computed, the cost as well as storage 
requirement for this k is doubled. For this reason, we would like to avoid negative j3 ik as much as possible. 
An extensive search performed in [29] gives the following preferred four stage, fourth order TVD Runge- 

Kutta method: 


u 




(4.15) 


u 


r n + 1 : 


louu zoiyouuu iww ■ ^ ^ 

53989 _ 102261 - 4806213 (1) 

2500000 5000000 1 ’ 20000000 


+ I tl ( 3) + lAtL(u< 3 >) 

3 o 

with a CFL coefficient c = 0.936 in (4.14). Notice that two L’s must be computed. The effective CFL 
coefficient, comparing with an ideal case without L’s, is 0.936 x | = 0.624. Since it is difficult to solve the 
global optimization problem, we do not claim that (4.15) is the optimal 4 stage, 4th order TVD Runge-Kutta 

method. 

A fifth order TVD Runge-Kutta method is also given in [69], 

For large scale scientific computing in three space dimensions, storage is usually a paramount considera- 
tion. There are therefore discussions about low storage Runge-Kutta methods [81], [10], which only require 
2 storage units per ODE equation. In [29], we considered the TVD properties among such low storage 
Runge-Kutta methods and found third order low storage TVD Runge-Kutta methods. 
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The general low-storage Runge-Kutta schemes can be written in the form [81], [10]: 

du {i) = Aidu {i -V + AfL(u ( ’ -1) ) 

( 4 -!6) « (<) =u^ + Bidu^\ i=l, ...,m 

tr (0) =u n , u^ = u n+1 , A 0 =0 

Only u and du must be stored, resulting in two storage units for each variable. 

Carpenter and Kennedy [10] have classified all the three stage, third order (m=3) low storage Runge- 
Kutta methods, obtaining the following one parameter family: 

z\ = + 36c| - 135c 2 + 84c 2 - 12 

z 2 = 2c^ + C2 — 2 

z 3 = 12 c\ - 18(4 + 18c| - llc 2 + 2 
z 4 = 36 c4 - 36(4 + - 8c 2 + 4 

z 5 = 69<4 - 62c^ + 28c 2 - 8 
z 6 = 34c^ - 46c| + 34c| - 13c 2 + 2 
(4.17) Bl = c 2 

g _ 12 c 2 (c 2 - l)(3z 2 — z \ ) - (3z 2 - zif 
2 144c 2 (3c2-2)(c 2 -1) 2 

B -24(3c 2 -2)(c 2 -1) 2 

(3 z 2 - z\) 2 - 12 c 2 (c 2 - 1)(3 z 2 - z\) 

A _ ~ z i(^ c 2 ~ 4c 2 + 1) -f 3 z 3 

2 ~~ (2 c 2 + 1)24 - 3(c 2 + 2)(2c 2 - l ) 2 

A = ~Z4Zi -f 1Q8(2 c 2 - 1 )cl - 3(2c 2 - 1 )z 5 

3 24z iC2 (c 2 - l ) 4 + 72c 2 z 6 + 72c|(2c 2 - 13) 

In [29] we converted this form into the form (4.8), by introducing three new parameters. Then we 
searched for values of these parameters that would maximize the CFL restriction, by a computer program. 
The result seems to indicate that 


(4.18) 


c 2 = 0.924574 


gives an almost best choice, with CFL coefficient c = 0.32 in (4.9). This is of course less optimal than (4.11) 
in terms of CFL coefficients, however the low storage form is useful for large scale calculations. 

We end this subsection by quoting the following numerical example [29], which shows that, even with a 
very nice second order TVD spatial discretization, if the time discretization is by a non-TVD but linearly 
stable Runge-Kutta method, the result may be oscillatory. Thus it would always be safer to use TVD 
Runge-Kutta methods for hyperbolic problems. 

The numerical example uses the standard minmod based MUSCL second order spatial discretization [79]. 
We will compare the results of a TVD versus a non-TVD second order Runge-Kutta time discretizations. 
The PDE is the simple Burgers equation 


(4.19) 

with a Riemann initial data: 



(4.20) 


u(x, 0) 


1, ifx<0 
—0.5, if x > 0 
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The nonlinear flux (\u 2 ) x in (4.19) is approximated by the conservative difference 

Ax (^ i+ * ' 

where the numerical flux fi+i is defined by 

f i+ , = h(u7 +v u+ +i ) 


with 


u7 = m + -minmod(ui+\ — Ui^Ui — Ui-i), 


i+ 


ut i = 1 — \minmod(ui+2 - Ui+i,Ui+i - U{) 

*+2 2 

The monotone flux h is the Godunov flux defined by (2.70), and the minmod function is given by 

minmod(a, b) = — mm(|a|, | 6 |). 

It is easy to prove, by using Harten’s Lemma [33], that the Euler forward time discretization with this second 
order MUSCL spatial operator is TVD under the CFL condition (4.6): 


(4.21) 


At < 


Ax 


2 ma Xj \u™\ 


Thus At = As t will be used in all our calculations. Actually, apart from a slight difference (the 

minmod function is replaced by a minimum-in-absolute-value function), this MUSCL scheme is the same as 

the second order ENO scheme discussed in Sect. 2.3.1. 

The TVD second order Runge-Kutta method we consider is the optimal one (4.10). The non- TVD 

method we use is: 


(4 22 ) u ^ = w n — 20A tL(u n ) 

u n+ 1 = u n + ^AtL(u") - ^AtL(u^). 

It is easy to verify that both methods are second order accurate in time. The second one (4.22) is however 
clearly non-TVD, since it has negative /3’s in both stages (i.e. it partially simulates backward in time with 

wrong upwinding). 

If the operator L is linear (for example the first order upwind scheme applied to a linear PDE), then 
both Runge-Kutta methods (actually all the two stage, second order Runge-Kutta methods) yield identical 
results (the two stage, second order Runge-Kutta method for a Unear ODE is unique). However, since our 
L is nonlinear, we may and do observe different results when the two Runge-Kutta methods are used. 

In Fig. 4.1 we show the result of the TVD Runge-Kutta method (4.10) and the non-TVD method (4.22), 
after the shock moves about 50 grids (400 time steps for the TVD method, 528 time steps for the non-TVD 
method). We can clearly see that the non-TVD result is osciUatory (there is an overshoot). 

Such oscillations are also observed when the non-TVD Runge-Kutta method coupled with a second order 
TVD MUSCL spatial discretization is apphed to a Unear PDE ( u t + u x = 0) (the scheme is still nonUnear 
due to the minmod functions). Moreover, for some Runge-Kutta methods, if one looks at the intermediate 
stages, i.e. u (i) for 1 < i < m in (4.8), one observes even bigger osciUations. Such oscillations may render 
difficulties when physical problems are solved, such as the appearance of negative density and pressure for 
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Fig. 4.1. Second order TVD MUSCL spatial discretization. Solution after 500 time steps. Left: TVD time discretization 
(f.10); Right: non-TVD time discretization (4-SS). 

Euler equations of gas dynamics. On the other hand, TVD Runge-Kutta method guarantees that each 
middle stage solution is also TVD. 

This simple numerical test convinces us that it is much safer to use a TVD Runge-Kutta method for 
solving hyperbolic problems. 

4.2.2. TVD multi-step methods. If one prefers multi-step methods rather than Runge-Kutta meth- 
ods, one can use the TVD high order multi-step methods developed in [66]. The philosophy is very similar 
to the TVD Runge-Kutta methods discussed in the previous subsection. One starts with a method of lines 
approximation (4.2) to the PDE (4.3), and an assumption that the first order Euler forward in time dis- 
cretization (4.4) is stable under a certain norm (4.5), with the time step restriction (4.6). One then looks 
for higher order in time multi-step methods such that the same stability result (4.5) holds, under a perhaps 
different restriction on At in (4.7), where c is again termed the CFL coefficient for the high order time 
discretization. 

The general form of the multi-step methods studied in [66] is: 


(4.23) 


.71+ 1 


= + At/3 k L{u n ~ k )) , 


Ac— 0 


Similar to the Runge-Kutta methods in the previous subsection, if all the coefficients are nonnegative 
a k > 0, /3k > 0, then (4.23) is just a convex combination of the Euler forward operators, with At replaced 
by since by consistency = 1. We thus have 

Lemma 4.3. [66] The multi-step method (4.23) is TVD under the CFL coefficient (4.7): 

(4.24) c=min^, 

* 0k 

provided that a* > 0, (3k > 0. 


In [66], schemes up to third order were found to satisfy the conditions in Lemma 4.3. Here we list a few 
examples. 

The following three step (m = 2) scheme is second order and TVD 

(4-25) u" +I = + \At L(u n ) + in"" 2 

4 2 4 

with a CFL coefficient c = 0.5 in (4.24). This translates to the same efficiency as the optimal second order 
TVD Runge-Kutta scheme (4.10), as here only one residue evaluation is needed per time step. Of course, 
the storage requirement is bigger here. There is also the problem of the starting values u 1 and u 2 . 
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The following five step (m = 4) scheme is third order and TVD 

(4.26) u n+l = + |A tL(u") + ^u"" 4 + ^AtL(u" -4 ) 

with a CFL coefficient c = 0.5 in (4.24). This translates to a better efficiency than the optimal third order 
TVD Runge-Kutta scheme (4.11), as here only one residue evaluation is needed per time step. Of course, 
the storage requirement is much bigger here. There is also the problem of the starting values u*, u 2 , u 3 and 
u 4 . 

There are many other TVD multi-step methods satisfying the conditions in Lemma 4.3 listed in [66]. It 
seems that if one uses more storage (larger m) one could get better CFL coefficients. 

In [66] we have been unable to find multi-step schemes of order four or higher satisfying the condition 
of Lemma 4.3. As in the Runge-Kutta case, we can relax the condition (3k > 0 by introducing the adjoint 
operator L . We thus have 

Lemma 4.4. [66] The multi-step method (4.23) is TVD under the CFL coefficient (4.7): 

(4.27) 

provided that a*, > 0, and L is replaced by L for negative /?*. 


Again, notice that, if we have both positive and negative /?*’ s, then both L(u n ) and L(u n ) must be 
computed, the cost as well as storage requirement will thus be doubled. 

We list here a six step (m = 5), fourth order multi-step method which is TVD with a CFL coefficient 
c = 0.245 in (4.24) [66]: 


,.n+l 


(4.28) 


1280 

256 


A tL{u n ) 


237 
4 + 


' 1) + io “ n ~ 5 


-AtL(u n ~ 5 ) 


4.2.3. The Lax-Wendroff procedure. Another way to discretize the time variable is by the Lax- 
Wendroff procedure [51]. This is also referred to as the Taylor series method for discretizing the ODE (4.2). 
We will again use the simple ID scalar conservation law (4.3) as an example to illustrate the procedure, 
however it applies to more general multidimensional systems. 

Starting from a Taylor series expansion in time: 

At 2 

(4.29) u(x,t + At) = u(x,t) + u t (x,t)At + u tt (x,t)— + ... 

The expansion is carried out to the desired order of accuracy in time. For example, a second order in time 
would need the three terms written out in (4.29). We then use the PDE (4.3) to replace the time derivatives 
by the spatial derivatives: 


U t(x,t) — /(u(x,f))x — / (ti(x, t)) U x {x, t), 

(4.30) u tt (x,t) = -(f(u(x,t)) tx = -{f'{u(x,t)u t (x,t)) x 

= (( f'{u{x,t)) 2 u x (x,t)) x 

= 2/ , (tt(x, t)) f"(u(x, t) (ttx(x, t)) 2 + (/ (u(x,t))) u xx (x,t); 
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This little exercise in (4.30) should convince us that it is always possible to write all the time derivatives 
as functions of the u(x,t) and its spatial derivatives. But the expression could be terribly complicated, 
especially for multidimensional systems. 

Once this is done, we substitute (4.30) into (4.29), and then discretize the spatial derivatives of u(x, £) 
by whatever methods we use. For example, in the cell averaged (finite volume) ENO schemes discussed 
in Sect. 2.3.1, we proceed as follows. We first integrate the PDE (2.65) in space-time over the region 
x [t n ,t n+1 ] to obtain 

(4.3i) ^ +1= ^-z h:{£ fwx i+i ,t))dt - £ 

Then, we use a suitable Gaussian quadrature to discretize the time integration for the flux in (4.31): 

( 4 - 32 ) At/n f(u(x i+ k,t))dt « ^2w a f{u(x i+ i,t n + 0 a At), 

a 

where (3 a and w a are Gaussian quadrature nodes and weights. Next we replace each 


f(u(x i+ i,t n + p a At) 


by a monotone flux: 

( 4 - 33 ) f(u(x i+ ^t n +/3 a At) &h(u(x~ + L,t n + /3 a At),u(x+ + ±,t n + /3 a At)') 

and use the Lax-Wendroff procedure (4.29)-(4.30) to convert 


U ( X t+l> tTl + A* At) 

to and its spatial derivatives also at t n y which can then be obtained by the reconstructions p(x) 

inside J* and J i+ i. Notice that the accuracy is just enough in this procedure, as each derivative of the 
reconstruction p{x) will be one order lower in accuracy, but this is compensated by the At in front of it in 
(4.29). 

This Lax-Wendroff procedure, comparing with the method of lines approach coupled with TVD Runge- 
Kutta or multi-step time discretizations, has the following advantages and disadvantages. 

Advantages: 

1. This is a truly one step method, hence it is quite compact (a second order method in space and time 
uses only three cells on time level n to advance to time level n -K 1 for one cell), and there are no 
complications such as boundary conditions needed in middle stages; 

2. It utilizes the PDE more extensively than the method of fines approach. This is also one reason that 
it can be so compact. 

Disadvantages: 

1. The algebra is very, very complicated for multi dimensional systems. This also increases operation 
counts for complicated nonlinear systems; 

2. It is more difficult to prove stability properties (e.g. TVD) for higher order methods in this frame- 
work; 

3. It is difficult and costly to apply this procedure to the conservative finite difference framework 
established in Sections 2.3 and 3.3. 
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4.3. Formulation of the ENO and WENO Schemes for the Hamilton- Jacobi Equations. In 

this section we describe high order ENO and WENO approximations to the Hamilton- Jacobi equation: 

f </>t + H(<f) X , (f>y ) = 0 

(4 - 34) { 4>(x,y,0) = 4>°(x,y) 

where H is a locally Lipschitz continuous Hamiltonian and the initial condition <p°(x, y) is locally Lipschitz 
continuous. We have written the equation (4.34) in two space dimensions, but the discussion is valid for 
other space dimensions as well. 

As is well known, solutions to (4.34) are Lipschitz continuous but may have discontinuous derivatives, 
regardless of the smoothness of <p°(x,y). The non-uniqueness of such generalized solutions also necessitates 
the definition of viscosity solutions, to single out a unique, practically relevant solution. The viscosity 
solution to (4.34) is a locally Lipschitz continuous function <j>(x, y,t), which satisfies the initial condition and 
the following property: for any smooth function ip(x,y,t ), if (x () , yo, to) is a local maximum point of <f> — v:, 
then 


(4.35) 


ipt{xo,yo,to) + H(ip x (xo,yo,to) + tM x o,2/o,<o)) ^ °> 


and, if (x 0 , yo,to) is a local minimum point of <j> - ip, then 


(4.36) 


i>t{x o,yo,to) + H(r{j x (xo,yo,to) + ^ y (x 0 ,2/o,*o)) > 0. 


Of course, the above definition means that whenever <f>(x, y, t ) is differentiable, (4.34) is satisfied in the 
classical sense. Viscosity solution defined this way exists and is unique. For details and equivalent definitions 
of viscosity solutions, see Crandall and Lions [22], 

H aTnil t.nn-Jar.ohi equations are actually easier to solve than conservation laws, because the solutions are 


typically continuous (only the derivatives are discontinuous). 

As before, given mesh sizes Ax, A y and At, we denote the mesh points as (x, , j/j . t n ) = ( iAx, j Ay , nAt). 
The numerical approximation to the viscosity solution 4>(xi,yj,t n ) of (4.34) at the mesh point ( Xi,yj,t n ) 


is denoted by <j>™y We again use a semi-discrete (discrete in the spatial variables only) formulation as a 
middle step in designing algorithms. In such cases, the numerical approximation to the viscosity solution 
<f>(xi,yj, t) of (4.34) at the mesh point (x,,yj, t) is denoted by the temporal variable t is not discretized. 

We will also use the notations D|<fe = and D y ± <j>ij = to denote the first order 

forward/ backward difference approximations to the left and right derivatives of < p(x , y) at the location (xi,yj). 

Since the viscosity solution to (4.34) is usually only Lipschitz continuous but not everywhere differen- 
tiable, the formal order of accuracy of a numerical scheme is again defined as that determined by the local 
truncation error in the smooth regions of the solution. Thus, a monotone scheme of the form 


(4.37) C/ 1 = G(tf-pj-ri ■ ’ ■,#+*,•+.) 

where G is a non-decreasing function of each argument, is called a first order scheme, although the provable 
order of accuracy in the norm is just \ [23], In the semi-discrete formulation, a five point monotone 
scheme (it does not pay to use more points for a monotone scheme because the order of accuracy of a 
monotone scheme is at most one [36]) is of the form 

(4.38) = ~H{D^(p l j(t), Dl^(t), £>* <M0, D y _<j> t} (t)). 
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The numerical Hamiltonian H is assumed to be locally Lipschitz continuous, consistent with H: H(u , u } v, v) = 
H(u,v), and is non-increasing in its first and third arguments and non-decreasing in the other two. Sym- 
bolically H(l, t, T)- ^ is easy to see that, if the time derivative in (4.38) is discretized by Euler forward 
differencing, the resulting fully discrete scheme, in the form of (4.37), will be monotone when At is suitably 
small. We have chosen the semi-discrete formulation (4.38) in order to apply suitable nonlinearly stable high 
order Runge-Kutta type time discretization, see Sect. 4.2. 

Semi-discrete or fully discrete monotone schemes (4.38) and (4.37) are both convergent towards the 
viscosity solution of (4.34) [23]. However, monotone schemes are at most first order accurate. As before, we 
will use the monotone schemes as building blocks for higher order ENO and WENO schemes. 

ENO schemes were adapted to the Hamilton- Jacobi equations (4.34) by Osher and Sethian [59] and 
Osher and Shu [60]. As we know now, the key feature of the ENO algorithm is an adaptive stencil high order 
interpolation which tries to avoid shocks or high gradient regions whenever possible. Since the Hamilton- 
Jacobi equation (4.34) is closely related to the conservation law (3.21), in fact in one space dimension 
they are exactly the same if one takes u = <f> x , it is not surprising that successful numerical schemes for 
the conservation laws (3.21), such as ENO and WENO, can be applied to the Hamilton- Jacobi equation 
(4.34). ENO and WENO schemes, when applied to Hamilton-Jacobi equations (4.34), can produce high 
order accuracy in the smooth regions of the solution, and sharp, non-oscillatory corners (discontinuities in 
derivatives) . 

There are many monotone Hamiltonians [23], [59], [60]. In this section we mainly discuss the following 

two: 

1. For the special case H(u ) t;) = f(u 2 ,v 2 ) where / is a monotone function of both arguments, such as 
the example H(u,v) = y/u 2 T n 2 , we can use the Osher-Sethian monotone Hamiltonian [59]: 

(4.39) H os (u + ,u-,v + } v-) = f{u 2 ,v 2 ) 

where, if / is a non-increasing function of u 2 , u 2 is implemented by 

(4*40) u 2 = (min(u~, 0)) 2 + (max(w + ,0)) 2 

and, if / is a non-decreasing function of u 2 , u 2 is implemented by 

(4.41) u 2 = (min(u+,0)) 2 -f (max(u’,0)) 2 

Similarly for v 2 . This Hamiltonian is purely upwind (i.e. when H(u,v) is monotone in u in the 
relevant domain [u”,u+] x [n~,u + ], only or u + is used in the numerical Hamiltonian according 
to the wind direction), and simple to program. Whenever applicable it should be used. This flux is 
similar to the Engquist- Osher monotone flux (2.71) for the conservation laws. 

2. For the general H we can always use the Godunov type Hamiltonian [5], [60]: 

(4.42) H (u + ,u ,n + ,n ) = e%tuei{u- ,u+) ex tvei{v~ ,v+) v ) 
where the extrema are defined by 


(4.43) 


e *^ix£/(a,6) 


{ m in a <u<b if a < b 

maxb< u<Ca if cl b 


Godunov Hamiltonian is obtained by attempting to solve the Riemann problem of the equation 
(4.34) exactly with piecewise linear initial condition determined by u* and v ± . It is in general not 
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unique, because in general min u max„ H{u, u) ^ max,., min u H(u, v) and interchanging the order of 
the two ext' s in (4,42) can produce a different monotone Hamiltonian. 

Godunov Hamiltonian is purely upwind and is the least dissipative among all monotone Hamiltonians 
[57], However, it might be extremely difficult to program, since in general analytical expressions for 
things like min u max, H{u, v) can be quite complicated. The readers will be convinced by doing the 
exercise of obtaining the analytical expression and programming H G for the ellipse in ellipse case in 
image processing where H(u,v) = Vau 2 + 2buv + cv 2 . For this case the Osher-Sethian Hamiltonian 
H os does not apply. 

We are now ready to discuss about higher order ENO or WENO schemes for (4.34). The framework is 
quite simple: we simply replace the first order scheme (4.38) by: 

(4.44) u-(t), v±(t), Vij (t)) 

where are high order approximations to the left and right x-derivatives of <j>{x, t/,t) at (a;*, 

(4.45) (t) = (if ,Vj,t) + 0(Ax r ) 

Similarly for vf (<). Notice that there is no cell-averaged version now. 

The key feature of ENO to avoid numerical oscillations is through the following interpolation procedure 
to obtain uf (t) and uf (t). These are just the same ENO procedure we discussed before in Sect. 2.2. We 
repeat it here with its own notations: 

ENO Interpolation Algorithm: Given point values f{xj), j - 0,±1,±2,--- of a (usually piecewise 
smooth) function /( x) at discrete nodes Xj, we associate an r-th degree polynomial Pf^ l/2 (x) with each 
interval [x^Xj+i], with the left-most point in the stencil as , constructed inductively as follows: 

(1) Pj+ 1/2 C-*- ) ~ f[ x j] P f[ X j ) Xj+l] (x — Xj), k min j, 

(2) If and Pj+y 2 { x ) are both defined, then let 

=/[x.(l-l),-”,X (i-l) ,] 

"min min 

b (l) = /[* fc ( , 7 l )_ 1 »-".x fc (i T ‘) + ,_ 1 ] 


(i) If |o®| > |6'|, then c® = 6 (0 and A® in = k^ in 1} - 1; otherwise c® = a® and fc® in = k^J 1 , 

(ii) p/+ 1/a (x) = p/fcu*) + c® nS-?‘ 1(x Xt) • 


In the above procedure /[•,••■,•] are the standard Newton divided differences, inductively defined as 
/[xi,x 2 , • • • ,x fc+1 ] = with /[xi] = /(x t ). 

ENO Interpolation Algorithm starts with a first degree polynomial Pfl\/ 2 ( x ) interpolating the function 
f(x) at the two grid points Xj and Xj+i. If we stop here, we would obtain the first order monotone scheme. 
When higher order is desired, we will in each step add just one point to the existing stencil, chosen from 
the two immediate neighbors by the size of the two relevant divided differences, which measures the local 
smoothness of the function / (x) . 

The approximations to the left and right x- derivatives of (f> are then taken as 


(4.46) 




- dx n± i/ 2tj 


( Xl y 
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where ^±i/2,j O 1 ) ‘ s obtained by the ENO Interpolation Algorithm in the x-direction, with y = yj and t both 
fixed, vfj are obtained in a similar fashion. The resulting ODE (4.44) is then discretized by an r-th order 
TVD Runge-Kutta time discretization in Sect. 4.2 to guarantee nonlinear stability. More specifically, the 
high order Runge-Kutta method we use in Sect. 4.2 will maintain TVD (total- variation- diminishing) or other 
stability properties, if these properties are valid for the simple first order Euler forward time discretization 
of the ODE (4.44). Notice that this is different from the usual linear stability requirement for the ODE 
solver. We thus obtain both nonlinear stability and high order accuracy in time. The second order (r = 2) 
and third order (r = 3) methods we use which has this stability property are given by (4.10) and (4.11), 
respectively. 

Time step restriction is taken as 


(4.47) 


At ( — max 
y Ax u,v 




l 

+ -r— max 

Ay u,v 



< 0.6 


where the maximum is taken over the relevant ranges of u, v. Here 0.6 is just a convenient number used in 
practice. This number should be chosen between 0.5 and 0.7 according to our numerical experience. 

WENO schemes can be used in a similar fashion for Hamilton- Jacobi equations [45], We will not present 
the details here. 


5. Applications. 

5.1. Applications to Compressible Gas Dynamics. One of the main application areas of ENO 
and WENO schemes is compressible gas dynamics. 

In 3D, Euler equations of gas dynamics are written as 


U t + f(U) x + g(U) y + h{U) z = 0 


where 


U = ( p,pu,pv,pw,E ), 


/ (U) = (pu, pu 2 + P, puv, puw, u{E + P)), 
g(U) = (pv, puv, pv 2 + P , pvw, v(E + P)), 
h(U) = ( pw , puw, pvw, pw 2 + P, w(E + P)). 

Here p is density, (u, v, w) is the velocity, E is the total energy, P is the pressure, related to the total energy 
E by 

E = + j-p(u 2 + v 2 + w 2 ) 

7-12 

with 7 = 1.4 for air. 

For the form of Navier-Stokes equations, for the eigenvalues and eigenvectors needed for the characteristic- 
wise ENO and WENO schemes, and for those equations appearing in curvilinear coordinates, see, e.g. [71]. 
We mention the following applications of ENO and ^VEN O schemes for compressible flow calculations: 
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1. Shock tube problem. This is a standard problem for testing codes for shock calculations. However, 
it is not the best test case for high order methods, as the solution structure is relatively simple 
(basically piecewise linear). The set-up is a Riemann type initial data: 


U(x, 0) 


U L if x < 0 
Ur if X > 0 


The two standard test cases are the Sod’s problem: 


(pL,qL,P L ) = (1,0,1); ( pR,qR,PR ) = (0.125,0,0.1) 


and the Lax’s problem: 

(, p L ,q L ,P L ) = (0.445,0.698,3.528); {p R ,q R ,P R ) = (0.5,0,0.571). 

We show the results of WENO (third order and fifth order) schemes for the Lax problem, in Fig. 
5.1. Notice that “PS” in the pictures means a way of treating the system cheaper than the local 
characteristic decompositions (for details, see [43]). “A” stands for Yang’s artificial compression [83] 
applied to these cases [43]. 



FIG. 5.1. Shock tube, Lax problem, density, (a): third order WENO; (b): fifth order WENO; (c): fifth order WENO with 
cheaper characteristic decomposition; (d): fifth order WENO with artificial compression. 

Wc can see from Fig. 5.1 that WENO perform reasonably well for these shock tube problem. The 
contact discontinuity is smeared more than the shock, as expected. Artificial compression helps 
sharpening contacts. For this problem, which is not the most demanding, the less expensive “PS” 
version of WENO work quite well. 

ENO schemes on this test case perform similarly. We will not give the pictures here. See [70]. 

2. Shock entropy wave interactions. This problem is very suitable for high order ENO and WENO 
schemes, because both shocks and complicated smooth flow feature co-exist. In this example, a 
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moving shock interacting with an entropy wave of small amplitude. On a domain [0,5], the initial 
condition is: 

p = 3.85714; u = 2.629369; P = 10.33333; when x < 0.5 

P — e -csln ( fcx ); u = 0; P = 1; when z > 0.5 

where e and k are the amplitude and wave number of the entropy wave, respectively. The mean 
flow is a pure right moving Mach 3 shock. If e is small compared to the shock strength, the shock 
will march to the right at approximately the non-perturbed shock speed and generate a sound wave 
which travels along with the flow behind the shock. At the same time, the perturbing entropy wave, 
after “going through” the shock, is compressed and amplified and travels approximately at the speed 
of u + c where u and c are the velocity and speed of the sound of the mean flow left to the shock. 
The amplification factor for the entropy wave can be obtained by linear analysis. 

Since the entropy wave here is set to be very weak relative to the shock, any numerical oscillation 
might pollute the generated waves (e.g. the sound waves) and the amplified entropy waves. In our 
tests, we take e = 0.01 and k — 13. The amplitude of the amplified entropy waves predicted by the 
linear analysis is 0.08690716 (shown in the following figures as horizontal solid lines). 

In Fig. 5.2, we show the result (entropy) when 12 waves have passed through the shock. It is clear 
that a lower order method (more dissipative) will damp the magnitude of the transmitted wave more 
seriously, especially when the waves are traveling more and more away from the shock. We can sec 
that, while fifth order WENO with 800 points already resolves the passing waves well, and with 
1200 points resolves the waves excellently, a second order TVD scheme (which is a good one among 
second order schemes) with 2000 points still shows excessive dissipation downstream. If we agree 
that fifth order WENO with 800 points behaves similarly as second order TVD with 2000 points, 
then there is a saving of a factor of 2.5 in grid points. This factor is per dimension , hence for a 3D 
time dependent problem the saving of the number of space-time grids will be a factor of 2.5 4 « 40, 
a significant saving even after factoring in the extra cost per grid point for the higher order WENO 
method. 

ENO schemes behave similarly for this problem. 

There is a two dimensional version of this problem, when the entropy wave can make an angle with 
the shock. The simulation results again show an advantage in using a higher order method, in Fig. 
5.3. Several curves are clustered in Fig. 5.3 around the exact solution, belonging to various fourth 
and fifth order ENO or WENO schemes. The circles correspond to a second order TVD scheme, 
which dissipates the amplitude of the transmitted entropy wave much more rapidly. 

3. Steady state calculations. This is important both in gas dynamics and in other fields of applications, 
such as in semiconductor device simulation, Sect. 5.3. For ENO or TVD schemes, the residue does 
not settle down to machine zero during the time evolution. It will decay first and then hang at the 
level of the local truncation errors. Presumably this is due to the fact that the numerical flux is 
not smooth enough (it is only Lipschitz continuous but not C 1 ). Although this is not satisfactory, 
it does not seem to affect the final solution (up to the truncation error level, which is how accurate 
the solution will be anyway). 

WENO schemes are much better in getting the residues to settle down to machine zeroes, due to 
the smoothness of their fluxes. 

In Fig. 5.4 we show the result of a one dimensional nozzle calculation. The residue in this case 
settles down nicely to machine zeros. Both fourth and fifth order WENO results are shown. 
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Fig. 5.2. ID shock entropy wave interaction. Entropy. Top: fifth order WENO with 800 points; middle: fifth order 
WENO with 1200 points; bottom: second order TVD with 2000 points. 

4. Forward facing step problem. This is a standard test case for high resolution schemes [82], However, 
second order methods usually already work well. High order methods might have some advantage 
in resolving the slip lines. 

The set up of the problem is the following: the wind tunnel is 1 length unit wide and 3 length units 
long. The step is 0.2 length units high and is located 0.6 length units from the left-hand end of the 
tunnel. The problem is initialized by a right-going Mach 3 flow. Reflective boundary conditions are 
applied along the walls of the tunnel and in-flow and out-flow boundary conditions are applied at 
the entrance (left-hand end) and the exit (right-hand end). For the treatment of the singularity at 
the corner of the step, we adopt the same technique used in [82], which is based on the assumption 
of a nearly steady flow in the region near the corner. 

In Fig. 5.5 we present the results of fifth order WENO and fourth order ENO with 242 x 79 grid 
points. 
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Fig. 5.3. 2D shock entropy wave interaction. Amplitude of amplified entropy waves. 800 points ( about 20 points per 
entropy wave length). 



Fig. 5.4. Density. Steady quasi-ID nozzle flow. 34 points. Left: fourth order WENO; right: fifth order WENO. 


5. Double Mach reflection. This is again a standard test case for high resolution schemes [82]. How- 
ever, second order methods usually again already work well. High order methods might have some 
advantage in resolving the flow below the Mach stem. 

The computational domain for this problem is chosen to be [0,4] x [0, 1], although only part of it, 
[0,3] x [0,1], is shown [82]. The reflecting wall lies at the bottom of the computational domain 
starting from x = Initially a right-moving Mach 10 shock is positioned at x = y = 0 and makes 
a 60° angle with the x-axis. For the bottom boundary, the exact post-shock condition is imposed 
for the part from x = 0 to x — ~ and a reflective boundary condition is used for the rest. At the 
top boundary of our computational domain, the flow values are set to describe the exact motion of 
the Mach 10 shock. See [82] for a detailed description of this problem. 

In Fig. 5.6 we present the results of fifth order WENO and fourth order ENO with 480 x 119 grid 
points. 

6. 2D shock vortex interactions. High order methods have some advantages in this case, as it resolves 
the vortex and the interaction better. 
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DENSITY WENO-LF-5 




Fig. 5.5. Flow post a forward facing step. Density: 242 X 79 grid points. Top: fifth order WENO; bottom, fourth order 
ENO. 



Fig. 5.6. Double Mach reflection. Density: 480 x 119 grid points. Top: fifth order WENO; bottom: fourth order ENO. 


58 






The model problem we use describes the interaction between a stationary shock and a vortex. The 
computational domain is taken to be [0,2] x [0,1]. A stationary Mach 1.1 shock is positioned at 
x = 0.5 and normal to the x-axis. Its left state is (p, u, v, P) = (1, 0,1). A small vortex is 

superposed to the flow left to the shock and centers at (x c ,y c ) = (0.25,0.5). We describe the vortex 
as a perturbation to the velocity (ix, v), temperature (T = £) and entropy (S = In of the mean 
flow and denote it by the tilde values: 


(5.1) 

u = ere a ( 1-r2) sin0 

(5.2) 

v = — ere ot ( 1 “ T ) cos 0 

(5.3) 

f _ (7 ~ l)e 2 e 2a ^ 1_T ^ 


4c*7 

(5.4) 

o 

II 


where r = — and r = \/( x — x c ) 2 -f (y — y c ) 2 - Here € indicates the strength of the vortex, a 
controls the decay rate of the vortex and r c is the critical radius for which the vortex has the 
maximum strength. In our tests, we choose e = 0.3, r c = 0.05 and oc — 0.204. The above defined 
vortex is a steady state solution to the 2D Euler equation. 

We use a grid of 251 x 100 which is uniform in y but refined in x around the shock. The upper and 
lower boundaries are intentionally set to be reflective. The results (pressure contours) are shown in 
Fig. 5.7 for a fifth order WENO with the cheap “PS” way of treating characteristic decomposition 
for the system. 



Fig. 5.7. 2D shock vortex interaction. Pressure. Fifth order WENO-LF-5-PS. 30 contours . (a) t=0.05. (b) t=0.20. (c) 
t=0.35. 

In [27], interaction of a shock with a longitudinal vortex is also investigated by the ENO method. 

7. How does the finite difference version of ENO and WENO handle non-rectangular domain? As we 
mentioned before, as long as the domain can be smoothly transformed to a rectangle, the schemes 
can be handily applied. 

We consider, as an example, the problem of a supersonic flow past a cylinder. In the physical space, 
a cylinder of unit radius is positioned at the origin on a x — y plane. The computational domain is 
chosen to be [0, 1] x [0, 1] on £-rj plane. The mapping between the computational domain and the 
physical domain is: 

(5-5) x = ( R x - ( R x - 1)£) cos(0(277 - 1)) 

(5.6) y = (it^ - (i^ - i)£) sin(0(277 - 1)) 

where we take R x = 3, Ry = 6 and 0 = Fifth order WENO and a uniform mesh of 60 x 80 in 
the computational domain are used. 
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The problem is initialized by a Mach 3 shock moving toward the cylinder from the left. Reflective 
boundary condition is imposed at the surface of the cylinder, i.e. £ = 1, inflow boundary condition 
is applied at £ = 0 and outflow boundary condition is applied at rj = 0, 1, 

We present an illustration of the mesh in the physical space (drawing every other grid line), and the 
pressure contour, in Fig. 5.8. Similar results are obtained by the ENO schemes but are not shown 

here. 


( a ) -4 -3 -2-10 1 (b) -4 -3 -2 -1 0 1 

Fig. 5.8. Flow past a cylinder, (a) Physical grid, (b) Pressure. WENO-LF-5. 20 contours 

8. Finally, we use the following problem to illustrate more clearly the power of high order methods. 
Consider the following idealized problem for the Euler equations in 2D: The mean flow is p = 1, P — 
1, and (u, v) = (1, 1) (diagonal flow). We add, to this mean flow, an isentropic vortex (perturbations 
in (u, v) and the temperature T — , no perturbation in the entropy S = ^): 

(6u,8v) = ^-e°- 5(1_r2) (-i/,x) 


PHYSICAL GRID 30x40 



PRESSURE 60x80 



8T = - 


(7 


877H 




ss = o, 


where ( x,y ) = (x - 5, y - 5), r 2 = x 2 + y 2 , and the vortex strength e = 5. 

Since the mean flow is in the diagonal direction, the vortex movement is not aligned with the mesh 
direction. 

The computational domain is taken as [0, 10] >t [0,10], extended periodically in both directions. This 
allows us to perform long time simulation without having to deal with a large domain. As we will 
see, the advantage of the high order methods are more obvious for long time simulations. 
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Fig. 5.9. Vortex evolution. Cut at x = 5. Solid: exact solution; circles : computed solution . Top: t = 50 (after 5 time 
periods ); bottom: t = 100 ( after 10 time periods). Left: second order TVD scheme; right: fifth order WENO scheme. 

It is clear that the exact solution of the Euler equation with the above initial and boundary conditions 
is just the passive convection of the vortex with the mean velocity. 

A grid of 80 2 points is used. The simulation is performed until t = 100 (10 periods in time). As can 
be seen from Fig. 5.9, fifth order WENO has a much better resolution than a second order TVD 
scheme, especially for the larger time t = 100. 

5.2. Applications to Incompressible Flows. In this section we consider numerically solving the 
incompressible Navier-Stokes or Euler equations 

U t + UU X + VUy = ft{U XX + Uyy) ~ X 

(5.7) V t + UV X + VVy = jl(v xx + Vyy) ~ Py 

U x + Vy = 0 

or their equivalent conservative form 

Ut T (u ) x “h (wv)j| = p(\t XX T Uyy) p x 

( 5 - 8 ) U t + (tit;)* + ( V 2 )y = p{V XX + Vyy) - Py 

U X -\~Vy — 0 

where (u,v) is the velocity vector, p is the pressure, \i > 0 for the Navier-Stokes equations and p = 0 for 
the Euler equations, using ENO and WENO schemes. We do not discuss the issue of boundary conditions 
here, thus the equation is defined on the box [0, 27r] x [0, 27t] with periodic boundary conditions in both 
directions. We choose two space dimensions for easy presentation, although our method is also applicable 
for three space dimensions. 

In some sense equations (5.7) are easier to solve numerically than their compressible counter- parts in 
Sect. 5.1, because the latter have solutions containing possible discontinuities (for example shocks and 
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contact discontinuities). However, the solution to (5.7), even if for most cases smooth mathematically, may 
evolve rather rapidly with time t and may easily become too complicated to be fully resolved on a feasible 
grid. Traditional linearly stable schemes, such as spectral methods and high-order central difference methods, 
are suitable for the cases where the solution can be fully resolved, but typically produce signs of instability 
such as oscillations when small scale features of the flow, such as shears and roll-ups, cannot be adequately 
resolved on the computational grid. Although in principle one can always overcome this difficulty by refining 
the grid, today’s computer capacity seriously restricts the largest possible grid size. 

As we know, the high resolution “shock capturing” schemes such as ENO and WENO are based on the 
philosophy of giving up fully resolving rapid transition regions or shocks, just to “capture” them in a stable 
and somehow globally correct fashion (e.g., with correct shock speed), but at the same time to require a 
high resolution for the smooth part of the solution. The success of such an approach for the conservation 
laws is documented by many examples in these lecture notes and the references. One example is the one 
and two dimensional shock interaction with vorticity or entropy waves [70], [71]. The shock is captured 
sharply and certain key quantities related to the interaction between the shock and the smooth part of the 
flow, such as the amplification and generation factors when a wave passes through a shock, are well resolved. 
Another example is the homogeneous turbulence for compressible Navier-Stokes equations studied in [71]. 
In one of the test cases, the spectral method can resolve all the scales using a 256 2 grid, while third order 
ENO with just 64 2 points can adequately resolve certain interesting quantities although it cannot resolve 
local quantities achieved inside the rapid transition region such as the minimum divergence. The conclusion 
seems to be that, when fully resolving the flow is either impossible or too costly, a “capturing” scheme such 
as ENO can be used on a coarse grid to obtain at least some partial information about the flow. 

We thus expect that, also for the incompressible flow, we can use high-order ENO or WENO schemes 
on a coarse grid, without fully resolving the flow, but still get back some useful information. 

A pioneer work in applying shock capturing compressible flow techniques to incompressible flow is by 
Bell, Colella and Glaz [6], in which they considered a second order Godunov type discretization, investigated 
the projection into divergence- free velocity fields for general boundary conditions, and discussed accuracy of 
time discretizations. Higher order ENO and WENO schemes for incompressible flows are extensions of such 
methods. 

We solve (5.8) in its equivalent projection form 



where P is the Hodge projection into divergence- free fields, i.e., if 


, then u x + v y = 0 and 


v y — u x = Vy — u x . See, e.g., [6]. For the current periodic case the additional condition to obtain a unique 
projection P is that the mean values of u and v are preserved, i.e., u(x, y)dxdy — J Q n f Q n u(x, y)dxdy 

311(1 C t € '( x ’ y) dxd v - 1 1 *(*> v) dxd v- 

We use N x and N y (even numbers) equally spaced grid points in x and y, respectively. The grid sizes 
are denoted by Ax = ^ and Ay = ^jf-, and the grid points are denoted by Xi — iAx and yj — jAy. The 
approximated numerical values of u and v at the grid point (xi,yj) are denoted by and Vij. 


We first describe the numerical implementation of the projection P. In the periodic case this is easily 
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achieved in the Fourier space. We first expand u and v using Fourier collocation: 


N V Ny JV, 

(5.10) u N (x,y) = ^ Ukie I{kx+lv) , v N (x,y) = ^3 

l = -?JLk=-2f 

where / = \/— T, lifci and v^i are the Fourier collocation coefficients which can be computed from the point 
values Uij and Vij, using either FFT or matrix- vector multiplications. The detail can be found in, e.g., [9], 
Derivatives, either by spectral method or by central differences, involve only multiplications by factors d% or 
df in (5.10) because e *( fcat +lv) are eigenfunctions of such derivative operators. For example, 


(5.11) 

for spectral derivatives; 

(5.12) 


JX 2/sin(^) 
ik ~ Ax 


d x k = Ik , df = II 

2/sin(^) 


Ay 


for the second order central differences which, when used twice, will produce the second order central 
difference approximation for w xx , and 


(5.13) 


dl = 




2 — cos(/cAx))(7 — cos(/cAx)) 


Ax 


2/y/q — cos(/Ay))(7 — cos(fAy)) 
Ay 


for the fourth order central differences which, when used twice, will produce the fourth order central difference 
approximation 3 ° Wl for w xx . High order filters, such as the exponential filter [55], 

[46]: 


(5.14) 


at = e~ a( *fe )2 


T v _ --“(it) 


2p 


where 2 p is the order of the filter and a is chosen so that e~ a is machine zero, can be used to enhance the 
stability while keeping at least 2p-th order of accuracy. This is especially helpful when the projection P is 
used for the under- resolved coarse grid with ENO methods. We use the fourth order projection (5.13) and 
the filter (5.14) with 2p = 8 in our calculations. This will guarantee third order accuracy (fourth order in 
L i) of the ENO scheme. We will denote this combination (the fourth order projection plus the eighth order 


filtering) by P4. To be precise, if 


-P 4 


and um and Vki are Fourier collocation coefficients of 


u and t?, then the Fourier collocation coefficients of u and v are given by 


(5.15) 


a k a i 


c^(dfu-dgf)) 

K ) 2 + «) 2 ’ 


J k°i 


~4fc(<ffu -rffcfi) 

(dfc) a + K) 2 


where <x£ and af are defined by (5.14) with 2 p = 8, and d% and ctf are defined by (5.13). 

Next we shall describe the ENO scheme for (5.8). Since (5.8) is equivalent to the non- conservative form 
(5.7), it is natural to implement upwinding by the signs of u and u, and to implement ENO equation by 
equation (the component version described in Sect. 2.3.5). The r-th order ENO approximation of, e.g., (u 2 ) x 
is thus carried out using the ENO Procedure 4.2. We mention a couple of facts needing attention: 

1. Take f(x) = u 2 (x, y) with y fixed. We start with the point values fi = f(xi ); 
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2. The stencil of the reconstruction is determined adaptively by upwinding and smoothness of f(x). It 
starts with either Xj or Xj+\ according to whether u > 0 or u < 0. 

There are two ways to handle the second derivative terms for the Navier-Stokes equations. One can 
absorb them into the convection part and treat them using ENO. For example, f(x) = u 2 (x,y) can be 
replaced by f(x) = u 2 {x,y) - fiu(x,y) x , where u(x,y) x itself can be obtained using either ENO or central 
difference of a suitable order. The remaining procedure for computing f(x) x would be the same as described 
above. Another simpler possibility is just to use standard central differences (of suitable order) to compute 
the double derivative terms. Our experience with compressible flow is that there is little difference between 
the two approaches, especially when the viscosity fx is small. 

In the above we have described the discretization for the spatial derivatives 



y = yj 


We then use the third order TVD (total variation diminishing) Runge-Kutta method (4.11) to discretize the 
resulting ODE: 


(5.17) 



— P 4 Lij 


obtaining: 


(5.18) 




\ (1) i 

) + 5 A < 

\ (2) 9 

) + 3 a <” 


Notice that wc have used the property P 4 o P 4 = P 4 in obtaining the discretization (5.18) from (5.17). 
This explicit time discretization is expected to be nonlinearly stable under the CFL condition 


(5.19) 


A t 


max 




Ax 2 


+ 


v)] s 1 


For small (i (which is the case we are interested in) this is not a serious restriction on At. 

We present some numerical examples in the following. 

Example 5.1: This example is used to check the third order accuracy of our ENO scheme for smooth 
solutions. We first take the initial condition as 


(5.20) u(x, y, 0) = - cos(x) sin(y), v(x, y, 0) = sin(x) cos (y) 
which was used in [6]. The exact solution for this case is known: 

(5.21) u{x,y,t) = - cos(x) sin(y)e _2,lt , v(x, y,t) = sin(x) cos(y)e 

We take Ax = A y - j; with N = 32,64, 128 and 256. The solution is computed up to t = 2 and the L 2 
error and numerical order of accuracy are listed in Table 5.1. For the fi = 0.05 case, we list results both 
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Table 5.1 

Accuracy of ENO Schemes for (12.2). 


m 

/* = 

0 

p = 0.05, central 

n = 0.05, ENO 

m 

Z /2 error 

order 

L 2 error 

order 

L 2 error 

order 

mi 



5.28(-4) 


4.87(-4) 


Ell 

G2ffl 


3.20(-5) 

4.04 

3.09 (-5) 

3.98 

Ell 

i 


1.93(-6) 

4.05 

1.89(-6) 

4.03 


2.28(-7) 


1.18 (-7) 

4.03 

1.16(-7) 

4.03 


N 

p = 0 

M = 

0.05, central 

^ = 

0.05, ENO 

Z/2 diff 

order 

error 

Z 2 diff 

order 

error 

L2 diff 

order 

error 

32 

1.14(-1) 






3.60(-2) 



64 

1.40(-2) 

3.02 

1.96(-3) 

2.78(-3) 

3.52 

2. 66 (-4) 

2.93(-3) 

3.62 


128 

1.46(-3) 

3.26 

1.69(-4) 

1.81 (-4) 

3.94 

1.26(-5) 

wmii 


1.18(-5) 

256 

l.ll(-4) 

3.77 

8.78(-6) 

1.09(-5) 


6.91 (-7) 

Hli 

| 

7.15(-7) 


with fourth order central approximation to the double derivative terms (central) and with ENO to handle 
the double derivative terms by absorbing them into the convection part (ENO). We can clearly observe fully 
third order accuracy (actually better in many cases because the spatial ENO is fourth order in the L\ sense) 
in this table. 

Example 5.2 : This is our test example to study resolution of ENO schemes when the grid is coarse. It 
is a double shear layer taken from [6]: 


(5.22) 


f tanh((y - 7r/2)/p) y < tt 
\ tanh((37r/2 - y)/p) y > n 


v(x,y, 0) = Jsin(x) 


where we take p = tt/ 15 and 5 = 0.05. The Euler equations (p = 0) are used for this example. The solution 
quickly develops into roll-ups with smaller and smaller scales, so on any fixed grid the full resolution is lost 
eventually. For example, the expensive run we performed using 512 2 points for the spectral collocation code 
(with a 18-th order filter (5.14)) is able to resolve the solution fully up to t = 8, Fig. 5.10, top left, as 
verified by the spectrum of the solution (not shown here), but begins to lose resolution as indicated by the 
wriggles in the vorticity contour at t = 10 (not shown here). On the other hand, the ENO runs with 64 2 
and 128 2 points produces smooth, stable results Fig. 5.10, top right and bottom left. In Fig. 5.10, bottom 
right, we show a cut at x = tt for v at t = 8. This gives a better feeling about the resolution in physical 
space. Apparently with these coarse grids the full structure of the roll-up is not resolved. However, when 
we compute the total circulation 


(5.23) cq = I uj(x,y)dxdy = / udx + vdy 

J n Jdn 

around the roll-up by taking SI = [§,^r] x [0,2tt] and using the rectangular rule (which is infinite order 
accurate for the periodic case) on the line integrals at the right-hand-side of (5.13), we can see that this 
number is resolved much better than the roll-up itself, Table 5.2. 

As an application of ENO scheme for incompressible flow, we consider the motion of an incompressible 
fluid, in two and three dimensions, in which the vorticity is concentrated on a lower dimensional set [31]. 
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Fig. 5.10. Double shear layer . Contours of vorticity. t = 8. Top left: spectral with 512 2 points; top right: ENO with 64 2 
points; bottom left: ENO with 128 2 points; bottom right: the cut at x = n of v , spectral method with 512 2 points, ENO method 
with 64 2 and with 128 2 points. 

Table 5.2 

Resolution of the Total Circulation. 


t 

2 

4 

6 

8 

10 

ENO 64 2 

0.87300 

3.07100 

7.16889 

9.88063 

10.90122 

ENO 128 2 

0.87452 

2.97810 

7.30999 

10.34414 

11.79418 

spectral 512 2 

0.87433 

2.98029 

7.28308 

10.46212 

11,85875 


Prominent examples are vortex sheets and vortex filaments in three dimensions, and vortex sheets, vortex 
dipole sheets and point vortices in two dimensions. 

In three dimensions, the equations are written in the form 

& + t>V£ - Vv £ = 0 

(5.24) Vxt/ = { 

V • v = 0 

where f(x,y, z,t) is the vorticity vector, and v(x,y, z,t) is the velocity vector. 

In a vortex sheet, £ is a singular measure concentrated on a two dimensional surface, while in a vortex 
filament, £ is a function concentrated on a tubular neighborhood of a curve. 

We use an Eulerian, fixed grid, approach, that works in general in two and three dimensions. In the 
particular case of the two dimensional vortex sheet problem in which the vorticity does not change sign, the 
approach yields a very simple and elegant formulation. 

The basic observation involves a variant of the level set method for capturing fronts, developed in [59]. 
The formulation we use here regularizes general ill-posed problems via the level set approach, using the 
idea that a simple closed curve which is the level set of a function cannot change its index, i.e. there is an 
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automatic topological regularization. This is very helpful for numerical calculations. The regularization is 
automatically accomplished through the use of dissipative schemes, which has the effect of adding a small 
curvature term (which vanishes as the grid size goes to zero) to the evolution of the interface. The formulation 
allows for topological changes, such as merging of surfaces. 

The main idea is to decompose f into a product of the form 

(5.25) £ = 

where P is a scalar function, typically an approximate 5 function. The variable y is a scalar function whose 
zero level set represents the points where vorticity concentrates, and 77 represents the vorticity strength 
vector. This decomposition is performed at time zero and is of course not unique. 

The observation is that once a decomposition is found, the following system of equations yields a solution 
to the Euler equations, replacing the original set of equations (5.24). 

(ft + = 0 

(5.26) rj t + vVr] — Vv 77 = 0 

Vxti = 

V • v = 0 


These equations have initial conditions 


< p ( 0 > •) = <po 

■) = ^0 


where </> 0 , 770 and P are chosen so that (5.25) holds at time t = 0. Notice that (5.25) and (5.26) imply 
that is orthogonal to 77 , and div(rj) — 0 . This is enforced in the initial condition and is maintained 
automatically by (5.25) and (5.26). 

When P is a distribution, such as a d function, approaching P with a sequence of smooth mollifiers P t 
yields a sequence of approximating solutions. This is the approach used in numerical calculations, since the 
5 function can only be represented approximately on a finite grid. The parameter e is usually chosen to be 
proportional to the mesh size. 

The advantage of this formulation, is that it replaces a possibly singular and unbounded vorticity function 
£, by bounded, smooth (at least uniformly Lipschitz) functions <p and 77 . Therefore, while it is not feasible 
to compute solutions of (5.24) directly, it is very easy to compute solutions of (5.26). 

In two dimensions, the vorticity is given by 

) 

\ u{t,x,y) f 

and hence the Euler equations are given by 


uJt + vVuj — 0 
curl(v) = uj 
div(v) = 0 


(5.27) 

(5.28) 
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Our formulation (5.26), becomes 


<Pt + vV<p = 0 

(5.29) V t + vVt) = 0 

curl(v) = P{<p) , r] 
div(v ) — 0 


where 7 ] is now a scalar. 

If the vortex sheet strength 7 ] does not change sign along the curve, it can be normalized to rj = 1 and 
the equations take on a particularly simple and elegant form: 

(5.30) <Pt + v(ip)V<p = 0 
where the velocity v((p) is given by 

(5.31) v = - ^ ^ A -1 P{<p) 

In this case, the vortex sheet strength along the curve is given by (see (5.33)). 

Example 5.3: Vortex Sheets in 2D. We consider the periodic vortex sheet in two dimensions, i.e. 

P{ip) = 8((p) in (5.31). The three dimensional case is defined in detail later. The evolution of the vortex sheet 
in the Lagrangian framework has been considered by various authors. Krasny [47], [48] has computed vortex 
sheet roll-up using vortex blobs and point vortices with filtering. Baker and Shelley [4] have approximated 
the vortex sheet by a layer of constant vorticity which they computed by Lagrangian methods. In the context 
of our approach, their approximation corresponds to approximating the 5 function by a step function. 

In our framework, we use a fixed Eulerian grid, and approximate (5.30) by the third order upwind ENO 
finite difference scheme with a third order TVD Runge-Kutta time stepping. At every time step, the velocity 
v is first obtained by solving the Poisson equation for the stream function \k: 

AS- = -P(<fi) 


with boundary conditions 

\k(x, ±1) = 0 

and periodic in x. This is done by using a second order elliptic solver FISHPAK. Once T is obtained, the 
velocity is recovered by v = (-$ y , $*) by using either ENO or central difference approximations (we do not 
observe major difference among the two: the results shown are those obtained by central difference). Once 
v is obtained, upwind biased ENO is easily applied to (5.30). 

The initial conditions are similar to the ones in [48], i.e given by a sinusoidal perturbation of a flat sheet: 

<Po(x>y) = y + 0.05 sin(Trx) 

The boundary condition for <p are periodic, of the form: 

<p{t,-l,y) = y>(M,y) 


-1) = v(t,x,l) -2 
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Fig. 5.11. Two dimensional vortex sheet simulation, t = 4. Top left: ENO with 128 2 points, 5 function width e = 12 Ax = 
Yq / Top right: ENO with 256 2 points, S function width t = 12Ax = Bottom left : ENO with 512 2 points, 5 function width 
e = 24Ax = ^ ; Bottom right: ENO with 1024 2 points, 6 function width t = 48Ax = 


The S function is approximated as in [61], [77] by 

(5.32) *«(*)=/ £( 1 + cos (^)) if ^ <e 

[ 0 otherwise 

For fixed e, there is convergence as Ax — > 0 to a smooth solution. One can then take e -+ 0. This two 
step limit is very costly to implement numerically. Our numerical results show that one can take e to be 
proportional to Ax, but convergence is difficult to establish theoretically. 

In Fig. 5.11, top left, we present the result at t = 4, of using ENO with 128 2 grid points with the 
parameter e in the approximate 5 function chosen as e = 12Ax. We use the graphic package TECPLOT to 
draw the level curve of tp = 0. Next, we keep c = 12 Ax but double the grid points in each direction to 256 2 , 
the result of t = 4 is shown in Fig. 5.11, top right. Comparing with Fig. 5.11, top left, we can see that 
there are more turns in the core at the same physical time when the grid size is reduced and the S function 
width € is kept proportional to Ax. One might wonder whether the core structure of Fig. 5.11, top right, 
is distorted by numerical error. To verify that this is not the case, we keep e = 12 x ^ ^ fixed , and 

reduce Ax, Fig. 5.11, bottom left and right. The three pictures overlay very well, the bottom two pictures 
in Fig. 5.11 are indistinguishable, indicating that the core structure is a resolved solution to the problem 
and convergence is obtained with fixed e. By reducing e for the more refined grids, more turns in the core 
can be obtained in shorter time (pictures not shown). 

The smoothing of the S function, and the third order truncation error in the advection step and the 
second order error in the inverse Laplacian are the only smoothing steps in our method. 

We now give the same example in three dimensions. We first sketch the algorithm for initializing and 
computing a periodic 3D vortex sheet, using (5.26). 

We let P(<p) = S((p) (in practice 5 is replaced by an approximation). The zero level set of ip is the vortex 
sheet T(s), parameterized by surface area s. The variable rjo is chosen to fit the initial vortex sheet strength. 
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For instance, given any smooth test function g 

it,g) = (voS{<po ),g) 

= J %(r 0 (s))3(r 0 (s))|^^|(is 

Thus, the initial vortex sheet strength is given by 

(5-33) 

To obtain the velocity vector, one introduces the vector potential A, where 

v — V x A, div(A) = 0 

and solves the Poisson equation 

(5.34) AA = -P{<p)y 

To ensure that div(A) = 0 , we require that div(r}) = 0 and that V<£ ■ 77 = 0 initially. It is easy to see that 
these equalities are maintained as t increases. 

The boundary conditions for the velocity are U 2 (x, ±1, z) = 0 and periodic in x and z. To obtain the 
boundary conditions for A = (A 1 ,A 2> A 3 ), we use the divergence free condition on A in addition to the 
velocity boundary condition. Thus, 

(5.35) A x {x,±l,z) = A 3 (x,±l,z) = 0 

d y A2{x, ±l,z) = 0 

and periodic in x,z. The Neumann condition requires the following compatibility condition 

J = 0 

Three dimensional runs are much more expensive than two dimensional runs, not only because the 
number of grid points increases, but also because there are now four evolution equations (for <p and 77 ), and 
three potential equations. We still use the third order ENO scheme coupled with the second order elliptic 
solver FISHPAK, with 64 3 grid points, and e is chosen as 6 Ax, which is the same in magnitude as that 
used in Fig. 5.11 of Example 5.3. The boundary conditions for (p are similar to the ones in two dimensions: 
periodic in all directions (module the linear term in y ). The vortex sheet strength vector 77 is periodic in all 
directions. 

We first verify whether we can recover the two dimensional results with the three dimensional setting. 
We use the initial condition 

<p o(x, y, z) = y + 0.05 sin( 7 rx) 

which is the same as that for Example 5.3, and choose a constant initial condition for r? as t)o(x, y, z) = (0, 0, 1). 
We observe exact agreement with our two dimensional results in Example 5.3, Fig. 5.11. Next, we consider 
the truly three dimensional problem with the initial condition chosen as 

ip 0 (x, y,z) = y + 0.05 sin( 7 rx) + 0.1 sin( 7 r.z) 

and r] is chosen as r*,(x,y, z) - (0, -O.Ittcos^), 1) which satisfies the divergence free condition as well as 
the condition to be orthogonal to V^. In Fig. 5.12, left, we show the level set of = 0 for t = 5. We can 
clearly see the roll up process and the three dimensional features. The cut at the constants z — 0 plane is 
shown in Fig. 5.12, right. 
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Fig. 5.12. Three dimensional vortex sheet simulation, t = 5. ENO with 64 3 points. S function width e = 6 Ax. Left: three 
dimensional level surface ; Right: z — 0 plane cut. 


5.3. Applications in Semiconductor Device Simulation. An interesting application area for ENO 
and WENO schemes is the equations in semiconductor device simulations. During the last decade, semi- 
conductor device modeling has attempted to incorporate general carrier heating, velocity overshoot, and 
various small device features into carrier simulation. The popular wisdom emerging from such concentrated 
study holds that global dependence of critical quantities, such as mobilities, on energy and/or temperature, 
is essential if such phenomena are to be modeled adequately. 

This gives rise to the various energy transport models, including the hydrodynamic model and the 
ET model, see, e.g. [41]. Unlike the earlier drift- diffusion models, which are basically parabolic, these 
new models contain significant transport effects [42], thus calling for discretization techniques suitable for 
hyperbolic problems. 

In this section we present two of such models. 

The first one is the hydrodynamic model. It is obtained by taking the first three moments of the 
Boltzmann equation. In the conservative format the hydrodynamic model is written as follows. Define the 
vector of dependent variables as 

(5.36) u = (n, < 7 , r, W ), 

where n is the electron concentration, p — (cr, r) is the momenta, and W is the total energy. The equations, 
in two dimensions, take the form 


(5.37) 

«t + h («)* + h{u) y = c(u) + G(u, <f>) + (0, 0, 0, V- (/cVT)) 

where 


(5.38) 

f r u \ — / a 2 , < 7 2 r 2 err 5<jW <r 2 + r 2 , 

m ’ 3 mn 2mn ’ mn ’ 3 mn ° 3 m 2 n 2 ' 

(5.39) 

t ( \ - l T aT 2 , T 2 , ur <T 2 N 5tW a 2 + T 2 

2U m' mn 3^mn 2 mn 3 mn T 3 m 2 n 2 

(5.40) 

, , , n ^ T W-W 0 . 

r ’ _ ’ r )> 

<P Ip lw 

(5.41) 

G(u) = (0, —enFi , —enFi, — enF • v). 


Here, F is the electric field, obtained by solving a Poisson’s equation: 

(5.42) F = 

(5.43) V- (eV</>) = -en - 77 , 4 . 
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Fig. 5.13, The one dimensional n"^-n-n~*" channel. Left: the doping ti^; Right: the velocity v, comparison of the HD 
model and the reduced HD model. 


where Tid is the doping (a given function which is typically discontinuous). 

The second model is the the energy transport model, written as 

(5.44) + f(u) x = g(u ) XX 4- h(u). 

In equation (5.44), 

nE 

(5.45) u={en,—), 

(5.46) /(«) = 4>'n (efx(E), n e (E ) + D{E)), 

(5.47) g(u) - ( nD{E ), nD E {E)), 

dE 

(5.48) h(u) = (0, en/i(£)(^) 2 + ~{n ~ n d )nD(E ) - n(— \ col i)). 

It can be shown that the left hand side defines a hyperbolic system, since the eigenvalues of f(u) are real, 
for all positive n and T. 

We first present one dimensional numerical results. The one dimensional n+-n-n + channel we simulate 
is a standard silicon diode with a length of 0.6^m, with a doping defined by n d = 5 x 10 17 cm -3 in [0,0.1] 
and in [0.5, 0.6], and n d = 2 x 10 15 cm~ 3 in [0.15,0.45], joined by smooth junctions (Fig. 5.13, left). The 
lattice temperature is taken as T 0 = 300 K. We apply a voltage bias of vbias = 1.5V. We use the full HD 
model; the relevant parameters can be found in [41]. In Fig. 5.13, right, we present the simulated velocity 
using the HD model. The dashed line shows the result computed with a reduced HD model by ignoring 
the transport effects. This type of reduced HD models are used quite often in engineering, as they tend to 
reduce the numerical difficulty when standard (not high resolution) schemes are used. However, we can see 
here that there is significant difference in the simulated results. 

We now present numerical simulation results for one carrier, two dimensional MESFET devices. The 
third order ENO shock-capturing algorithm with Lax-Friedrichs building blocks, as described elsewhere in 
these lecture notes, is applied to the hyperbolic part (the left hand side) of Equations (5.37) and (5.44). 
The TVD third order Runge-Kutta time discretization (4.11) is used for the time evolution towards steady 
states. The forcing terms on the right hand side of (5.37) and (5.44) are treated in a time consistent way in 
the Runge-Kutta time stepping. The double derivative terms on the right hand side of (5.37) and (5.44) are 
approximated by standard central differences owing to their dissipative nature. The Poisson equation (5.43) 
is solved by direct Gauss elimination for one spatial dimension and by Successive Over-Relaxation (SOR) 
or the Conjugate Gradient (CG) method for two spatial dimensions. Initial conditions are chosen as n = n d 
for the concentration, T = T 0 for the temperature, and u = v = 0 for the velocities. A continuation method 
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is used to reach the steady state: the voltage bias is taken initially as zero and is gradually increased to the 
required value, with the steady state solution of a lower biased case used as the initial condition for a higher 
one. 

We simulate a two dimensional MESFET of the size 0.6 x 0.2 pm 2 . The source and the drain each 
occupies 0.1pm at the upper left and the upper right, respectively, with the gate occupying 0.2pm at the 
upper middle (Fig. 5.14, top left). The doping is defined by rid = 3 x 10 17 cm -3 in [0,0.1] x [0.15,0.2] and 
in [0.5, 0.6] x [0.15,0.2], and rid = lx 10 17 cm -3 elsewhere, with abrupt junctions (Fig. 5.14, top right). 
A uniform grid of 96 x 32 points is used. Notice that even if we may not have shocks in the solution, the 
initial condition n — rid is discontinuous, and the final steady state solution has a sharp transition around 
the junction. With the relatively coarse grid we use, the non- oscillatory shock capturing feature of the ENO 
algorithm is essential for the stability of the numerical procedure. 

We apply, at the source and drain, a voltage bias vbias — 2V. The gate is a Schottky contact, with a 
negative voltage bias vgate = — 0.8F and a very low concentration value n — 3.9 x 10 5 cm -3 . The lattice 
temperature is taken as To = 300° K. The numerical boundary conditions are summarized as follows (where 
= ^r~ l n kb = 0.138 x 10 -4 , e = 0.1602, and n* = 1.4 x 10 lo cm -3 in our units): 

• At the source (0 < x < 0.1 , y = 0.2): $ = $o for the potential; n — 3 x 10 17 cm“ 3 for the concen- 

tration; T = 300°AT for the temperature; u = Opm/ps for the horizontal velocity; and Neumann 
boundary condition for the vertical velocity v (i.c. = 0 where ft is the normal direction of the 

boundary) . 

• At the drain (0.5 < x < 0.6, y = 0.2): 3> = 3>o + vbias = 4>o + 2 for the potential; n =■ 3 x 10 17 cm~ 3 
for the concentration; T — 300 ° K for the temperature; u ~ 0 pm/ps for the horizontal velocity; and 
Neumann boundary condition for the vertical velocity v. 

• At the gate (0.2 < x < 0.4, y — 0.2): = $o + vgate = 3>o~0.8 for the potential; n = 3.9 x 10 5 cra~ 3 

for the concentration; T = 300 °K for the temperature; u = 0 pm/ps for the horizontal velocity; and 
Neumann boundary condition for the vertical velocity v. 

• At all other parts of the boundary (0.1 < x < 0.2, y = 0.2; 0.4 < x < 0.5, y = 0.2; x = 0, 0 < y < 0.2; 
x — 0.6, 0 < y < 0.2; and 0 < x < 0.6, y = 0), all variables are equipped with Neumann boundary 
conditions. 

The boundary conditions chosen are based upon physical and numerical considerations. They may not be 
adequate mathematically, as is evident from some serious boundary layers observable in the concentration (see 
pictures in [41]). ENO methods, owing to their upwind nature, are robust to different boundary conditions 
(including over-specified boundary conditions) and do not exhibit numerical difficulties in the presence of 
such boundary layers, even with the extremely low concentration prescribed at the gate (around 10“ 12 
relative to the high doping). We point out, however, that boundary conditions affect the global solution 
significantly. We have also simulated the same problem with different boundary conditions, for example with 
Dirichlet boundary conditions everywhere for the temperature, or with Neumann boundary conditions for 
all variables except for the potential at the contacts. The numerical results (not shown here) are noticeably 
different. This indicates the importance of studying adequate boundary conditions, from both a physical 
and a mathematical point of view. 

The velocity vectors resulting from the hydrodynamic model simulation are presented in Fig. 5.14, 
bottom left. In Fig. 5.14, bottom right, we compare the temperature at y = 0.175 from the simulations of 
the hydrodynamic model and of the ET model. Clearly there is a significant difference between these two 
models for this 2D case. 
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Fig. 5.14. Two dimensional MESFET. Top left: the geometry; Top right: the doping n<i; Bottom left: the velocity (u,v) 
obtained by the HD model; Bottom right: comparison of the hydrodynamic (HD) model (solid line) and the energy transport 
(ET) model (plus symbols), cut at the middle of the high doping blobs y = 0.175, the temperature T. 


There are also new models in semiconductor device simulation (e.g. [14], [15]), which are worthy of 
investigations. ENO and WENO schemes provide robust and reliable tools for carrying out such investiga- 
tions. 
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